S3WU> 

fIS’i 


NASA Contractor Report CR- 1871 06 

Investigation of Advanced Counterrotation 
Blade Configuration Concepts for High Speed 
Turboprop Systems 

Task II - Unsteady Ducted Propfan Analysis 
Final Report 


Edward J. Hall, Robert A. Delaney, and Janies L. Bettner 
Allison Gas Turbine Division of General Motors 
Indianapolis , Indiana 


May 1991 



Prepared for 

Lewis Research Center 

Under Contract NAS3-25270 


(NASA-CR-187106) INVESTIGATION OF ADVANCED 
COUNTERROTATION BLADE CONFIGURATION CONCEPTS 
FOR HIGH SPEED TURBOPROP SYSTEMS. TASK 2 ' 
UNSTEADY DUCTFD PROPFAN ANALYSIS Final 
Report (General Motors Corp.) 132 p 


N91-251 13 


Unc 1 as 

G 3/02 0020327 


PRECEDING PAGE BLANK NOT FILMED 


^_JLLjN*N*IO«AI» mm iii 



Preface 


This manual was prepared by Edward J. Hall, Robert A. Delaney, and James 
L. Bettner of the Allison Gas Turbine Division, General Motors Corporation, In- 
dianapolis, IN. The work was performed under NASA Contract NAS3-25270 from 
March, 1990 to March, 1991. The grid generation, flow code theory, and program- 
ming modifications necessary for the analysis of ducted propfans were performed by 
Edward J. Hall. The Allison program manager for this contract was James L. Bettner. 
The NASA program manager for this contract was Christopher J. Miller. 


Acknowledgements 

The authors would like to express their appreciation to the following NASA 
personnel who contributed to this program: 

Mark Celestina for his helpful suggestions concerning the development of the 
viscous solver. 

Dr. John J. Adamczyk for the many helpful technical discussions concerning the 
development of the computer codes. 

Dr. Cristopher J. Miller for his suggestions and critical review of the program. 

Dr. David P. Miller for the use of his grid generation code, TIGGERC 

Christopher E. Hughes for his helpful comments and suggestions 

The services of. the NASA Numerical Aerodynamic Simulation (NAS) facilities 





IV 


and personnel are gratefully acknowledged. 

UNIX is a trademark of AT&T 
IRIS is a trademark of Silicon Graphics, Inc. 
UNICOS is a trademark of Cray Computers 
PostScript is a trademark of Adobe Systems, Inc. 



TABLE OF CONTENTS 


NOTATION xv 

1. SUMMARY i 

2. INTRODUCTION 3 

3. GRID GENERATION ALGORITHM 9 

3.1 Computational Domain 9 

3.2 Unducted Propfan Grid Generation 12 

3.3 Ducted Propfan Grid Generation 18 

4. 3D EULER/NAVTER-STOKES NUMERICAL ALGORITHM . 29 

4.1 Nondimensionalization 29 

4.2 Governing Equations 30 

4.3 Runge-Kutta Time Integration 34 

4.4 Fluid Properties 37 

4.5 Turbulence Model 37 

4.6 Artificial Dissipation 41 

4.7 Implicit Residual Smoothing 43 

4.8 Boundary Conditions 46 



vi 

4.9 Multiple- Block Coupling 4g 

4.10 Solution Procedure 49 

5. RESULTS 51 

5.1 SR7 2-Bladed Propfan Modane Tests 52 

5.2 SR7 8-Bladed Propfan gg 

5.3 Ducted Propfan Test Case gQ 

5.4 NASA 1.15 Pressure Ratio Fan g7 

6. CONCLUSIONS 109 

REFERENCES m 

APPENDIX A. ADPAC DISTRIBUTION LIST 115 



vii 


LIST OF FIGURES 

Figure 2.1: Ducted propfan aerodynamic characteristics 6 

Figure 3.1: Ducted propfan analysis computational domain 11 

Figure 3.2: Meridional plane projection grid generation subregions for an 

unducted propeller 13 

Figure 3.3: Sample grid for an unducted propfan 19 

Figure 3.4: Axisymmetric plane projection grid generation subregions for 

a ducted propeller 21 

Figure 3.5: Sample multiple-block C-grid for a ducted propfan 28 

Figure 4.1: Three-dimensional finite volume cell 35 

Figure 5.1: SR7 propfan design characteristics 53 

Figure 5.2: 2-bladed SR7 propfan geometry and steady state inviscid cal- 
culation grid for Modane test comparison 56 

Figure 5.3: 2-bladed SR7 propfan geometry and steady state viscous cal- 
culation grid for Modane test comparison 57 

Figure 5.4: Comparison of predicted and experimental airfoil surface static 

pressure coefficient distributions for 2-bladed SR7 propfan Modane 
test (28.4% span) 58 



Vlll 


Figure 5.5: Comparison of predicted and experimental airfoil surface static 

pressure coefficient distributions for 2-bladed SR7 propfan Modane 
test (93.1% span) 59 

Figure 5.6: Comparison of predicted and experimental airfoil suction sur- 

face static/total pressure ratio contours for 2-bladed SR7 prop- 
fan 60 

Figure 5.7: Illustration of predicted turbulent flow leading edge vortex 

particle trajectory traces for for 2-bladed SR7 propfan ... 61 

Figure 5.8: Comparison of laminar and turbulent predicted suction surface 

shear flow patterns and static/total pressure ratio contours for 
2-bladed SR7 propfan 62 

Figure 5.9: Comparison of predicted and experimental airfoil surface static/total 

pressure ratio distributions for 2-bladed SR7 propfan (28.4% 
span, M=0.5, J=3.06) 65 

Figure 5.10: Comparison of predicted and experimental airfoil surface static/total 
pressure ratio distributions for 2-bladed SR7 propfan (94.4% 
span, M=0.5, J=3.06) 66 

Figure 5.11: Comparison of predicted and experimental airfoil surface time- 
dependent static/total pressure ratio history for 2-bladed SR7 
propfan (suction side, 4.9% chord, 64.0% span, 3 degrees angle 
of attack) 69 



Figure 5.12: Comparison of predicted and experimental airfoil surface time- 
dependent static/total pressure ratio history for 2-bladed SR7 
propfan (suction side, 36.7% chord, 64.0% span, 3 degrees an- 
gle of attack) 

Figure 5.13: Comparison of predicted and experimental airfoil surface time- 
dependent static/total pressure ratio history for 2-bladed SR7 
propfan (pressure side, 4.9% chord, 64.0% span, 3 degrees an- 
gle of attack) 

Figure 5.14: Comparison of predicted and experimental airfoil surface time- 
dependent static/total pressure ratio history for 2-bladed SR7 
propfan (pressure side, 10.0% chord, 64.0% span, 3 degrees 
angle of attack) 

Figure 5.15: Comparison of predicted and experimental airfoil surface time- 
dependent static/total pressure ratio history for 2-bladed SR7 
propfan (pressure side, 63.3% chord, 64.0% span, 3 degrees 
angle of attack) 

Figure 5.16: Comparison of predicted and experimental airfoil surface time- 
dependent static/total pressure ratio history for 2-bladed SR7 
propfan (suction side, 27.9% chord, 91.0% span, 3 degrees an- 
gle of attack) 



X 


Figure 5.17: Comparison of predicted and experimental airfoil surface time- 
dependent static/total pressure ratio history for 2-bladed SR 7 
propfan (suction side, 69.8% chord, 91.0% span, 3 degrees an- 
gle of attack) 75 

Figure 5.18: Comparison of predicted and experimental airfoil surface time- 
dependent static/total pressure ratio history for 2-bladed SR 7 
propfan (pressure side, 27.9% chord, 91.0% span, 3 degrees 
angle of attack) yg 

Figure 5.19: Comparison of predicted and experimental airfoil surface time- 
dependent static/total pressure ratio history for 2-bladed SR 7 
propfan (pressure side, 89.8% chord, 91.0% span, 3 degrees 
angle of attack) 77 

Figure 5.20: Instantaneous propfan surface static/total pressure ratio con- 
tours and blade tip particle trajectories for 2-bladed SR 7 prop- 
fan at angle of attack (M=0.5, J=3.06, angle of attack = 3 
degrees) 7 g 

Figure 5.21: Full rotor grid for SR7 time-dependent inviscid calculations . 81 

Figure 5.22. Predicted instantaneous propfan inviscid blade passage static 
pressure contours for 8 -bladed SR7 propfan at angle of attack 
(M=0.8) g 2 

Figure 5.23: Predicted instantaneous propfan surface static/total pressure 
ratio contours for 8 -bladed SR7 propfan at angle of attack 
(M= 0 . 8 , angle of attack = 4.6 degrees) g 3 



XI 


Figure 5.24: Predicted rotational single blade power coefficient histories for 
8-bladed SR7 propfan at angle of attack (M=0.8, angle of at- 
tack = 4.6 degrees) 84 

Figure 5.25: Full rotor grid system for ducted SR7 propfan geometry ... 88 

Figure 5.26: Instantaneous blade passage static pressure contours for ducted 

SR7 propfan geometry (M=0.8) 89 

Figure 5.27: Comparison of time-average and unsteady blade section static 

pressure distributions for ducted SR7 propfan (90% span). . . 90 

Figure 5.28: Predicted surface static pressure contours for ducted SR7 prop- 

fan geometry at angle of attack (M=0.8) 91 

Figure 5.29: Predicted single blade power coefficient rotational histories for 
ducted and unducted 8-bladed SR7 propfans at angle of attack 


(M=0.8) 92 

Figure 5.30: NASA 1.15 pressure ratio fan stage geometry (dimensions in 

cm) 94 


Figure 5.31: NASA 1.15 pressure ratio fan viscous flow grid system .... 95 

Figure 5.32: Predicted cowl surface velocity vectors in the vicinity of the 
outer surface shock for the NASA 1.15 pressure ratio fan stage 
geometry (M=0.85) 96 

Figure 5.33: Comparison of viscous and inviscid predicted and experimen- 
tal cowl surface leading edge static pressure ratio coefficient 
distributions for NASA 1.15 pressure ratio fan (M=0.75) . . 98 



Xll 


Figure 5.34: Comparison of viscous and inviscid predicted and experimen- 
tal cowl surface leading edge static pressure ratio coefficient 
distributions for NASA 1.15 pressure ratio fan (M=0.85) . . 99 

Figure 5.35: Full rotor grid system for NASA 1.15 pressure ratio fan . 100 

Figure 5.36: Predicted viscous instantaneous static pressure contours for 

NASA 1.15 pressure ratio fan at angle of attack. (M=0.20) . 101 

Figure 5.37: Comparison of viscous predicted time-averaged and unsteady 
envelope blade surface static pressure ratio distributions for 
NASA 1.15 pressure ratio fan at angle of attack (10% span, 
M=0 - 20 ) 102 

Figure 5.38: Comparison of viscous predicted time-averaged and unsteady 
envelope blade surface static pressure ratio distributions for 
NASA 1.15 pressure ratio fan at angle of attack (50% span, 

M=0 - 20 ) 103 

Figure 5.39: Comparison of viscous predicted time-averaged and unsteady 
envelope blade surface static pressure ratio distributions for 
NASA 1.15 pressure ratio fan at angle of attack (90% span, 
M=0.20) 2Q4 

Figure 5.40: Comparison of experimental and viscous predicted instanta- 
neous fan face total pressure ratio contours for NASA 1.15 
pressure ratio fan (M=0.2, angle of attack = 40 degrees). . . 106 



Figure 5.41: Comparison of experimental and instantaneous viscous pre- 
dicted cowl windward side internal static/inlet total pressure 
ratio distribution for NASA 1.15 pressure ratio fan (M=0.2, 
angle of attack = 40 degrees) 



xiv 


LIST OP TABLES 


Table 5.1: 


Summary of computational characteristics for ADPAC test cases 52 



XV 


NOTATION 


A list of the symbols used throughout this document and their definitions is 
provided below for convenience. Parameter values are shown in parentheses. 

Roman Symbols 


a . . . speed of sound 

cp . . . specific heat at constant pressure 

Cy . . . specific heat at constant volume 

e . . . internal energy 

i . . . z index of numerical solution 

j ... r index of numerical solution 

k . . . 0 index of numerical solution or thermal conductivity 
/ . . . turbulence model damping function 

n... time step index of numerical solution or rotational speed (revolutions/sec) 
n . . . outward unit normal vector 
p . . . pressure 

r . . . radius or radial coordinate 
t . . . time 
v... velocity 
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z ... axial coordinate 

A. . . surface area 

-A'*’ . . . turbulence model parameter (26) 

ADPAC . . . Advanced Ducted Propfan Analysis Codes 
AO A . . . Angle of Attack aerodynamic analysis code 
AOAPLOT . . . Ducted propfan automated plotting program 
ASCII . . . American Standard Code for Information Interchange 

B . .. number of propeller blades 

Cp... power coefficient (Cp = P/pn^D^) 

Ci-.. thrust coefficient ( Ci = T/pn^D^) 

Cep . . . turbulence model parameter (1.6) 

C kleb • * * turbulence model parameter (0.3) 

C wake • • • turbulence model parameter (0.25) 

CFL ... Courant-Freidrichs-Levy number (At/At max 5 ^ a j/ e ) 

CHGRIDV 2 . . . Ducted propfan grid generation code 

D . . . dissipation flux vector, turbulent damping parameter, or diameter 

DELT . . . adjacent cell grid spacing 

DELTLE . . . adjacent cell grid spacing at leading edge 

DELTT E ... adjacent cell grid spacing at trailing edge 

F . . . flux vector in z direction or turbulence model function 

G . . . flux vector in r direction 

H . .. flux vector in 6 direction 
Ht ... total enthalpy 
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J . .. advance ratio (J = U/nD) 

K . .. source term flux vector or turbulence model parameter (0.0168) 

L . . . length 

M . . . Mach number 

O . .. orthogonality 

P . .. power 

Pr . . . Prandtl number 

turbulent • • • turbulent Prandtl number (0.9) 

Q . . . vector of dependent variables 

R. . . gas constant or residual or maximum radius 
RAT . . . maximum ratio of adjacent cell grid spacings 
ROTCGRID . . . Ducted propfan full rotor grid rotation program 
ROTCF LOW . . . Ducted propfan full rotor flow rotation program 
SDBLI B . . . Scientific DataBase Library (binary file formats) 

5 ... arc length or pertaining to surface area normal 

T . .. temperature or torque 

U . . . flight velocity 

V . . . volume 


Greek Symbols 


a. .. time-stepping factor 

P ... local propfan blade angle 

^3 ... 3/4 radius propfan blade setting angle 
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c . . . modified second-order damping coefficient 

... modified fourth-order damping coefficient 

p ... density 
o 

k ... second-order damping coefficient 
. . . fourth-order damping coefficient 
7 . . . specific heat ratio 

6 .. . spatial second-order central difference operator 
X... blockage factor 

X v .. . second coefficient of viscosity (= — j/t) 
p ... coefficient of viscosity 
Tf ... radial transformed variable 
£ ... axial transformed variable 
C • • • circumferential transformed variable 
v ... damping factor 
T . . . boundary layer dissipation factor 
A... increment of change 

Special Symbols 

V . . . spatial vector gradient operator 
A . . . spatial forward difference operator 
V • • • spatial backward difference operator 


Superscripts 
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[ ] . . . averaged variable 

[ ] . • . dimensional variable 

[ ] . . . implicitly smoothed variable 

[ ] . . . vector variable 

[ ]* . . . intermediate variable 

[ ] n . . . time step index of variable 

Subscripts 

[ \ e j f ec ti ve • • • effective flow value 
[ ]ij y k • • • grid point index of variable 
[ haminar • • • laminar flow value 
[ }max • • • maximum value 
[ } m i n . . . minimum value 
[ ]p . . . related to pressure 
[ ]ps • • • pressure (high pressure) surface 
[ ]ss • • • suction (low pressure) surface 
[ ]t ■ • • total quantity 

[ ]z • • • derivative or value with respect to z 
[ ]r • • • derivative or value with respect to r 
[]$... derivative or value with respect to 9 
[ ] turbulent * * * turbulent flow value 
[ ]oo • • • freestream value 
[ ] re f . . . reference value 
[ ]&/ e j • • • Klebaiioff intermit ten cy factor 
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[ ] W ake • • * turbulent flow wake parameter 
[ ]2 • • • second-order value 
[ ]4 . . . fourth-order value 



1. SUMMARY 


The primary objective of this study was the development of a time-dependent 
three-dimensional Euler/Navier-Stokes aerodynamic analysis to predict unsteady com- 
pressible transonic flows about ducted and unducted propfan propulsion systems at 
angle of attack. The computer codes resulting from this study are part of a group of 
codes referred to as ADPAC (Advanced Ducted Propfan Analysis Codes). This doc- 
ument is the final report detailing the development and application of the ADPAC 
codes developed under Task II of NASA Contract NAS3-25270, Unsteady Ducted 
Propfan Analysis. 

Aerodynamic calculations were based on a four-stage Runge-Kutta time-marching 
finite volume solution technique with added numerical dissipation. A time-accurate 
implicit residual smoothing operator was utilized for unsteady flow predictions. A 
single H-type grid was used to discretize each blade passage for unducted propfans. 
For ducted propfans, a coupled system of five grid blocks utilizing an embedded C- 
grid about the cowl leading edge were used to discretize each blade passage. Grid 
systems were generated by a combined algebraic/elliptic algorithm developed specif- 
ically for ducted propfans. Numerical calculations were compared with experimental 
data for both ducted and unducted propfan flows. The solution scheme demonstrated 
efficiency and accuracy comparable with other schemes of this class. 
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2. INTRODUCTION 


The development and demonstration of propfan propulsion systems in the last 
decade can be cited as a prime example of aerodynamic technology which has directly 
benefitted from the use of computational fluid dynamics (CFD) tools (see Hager 
and Vrabel [1]). The highly three-dimensional nature of the blade designs do not 
lend themselves easily to standard two-dimensional design philosophies, and therefore 
require advanced CFD analysis techniques to verify the performance of candidate 
designs. Numerical analyses for steady propfan flows in three spatial dimensions have 
been given by Barton et al. [2], Usab et al. [3], Saito et al. [4], Matsuo et al. [5], 
Kobayawa and Hatano [6], Celestina et al. [7], Whitfield et al. [8] and others. 

Recently, the concept of utilizing a ducted propfan (ultra high bypass fan) as a 
primary propulsion system has been considered, and has also been an area of intense 
CFD study. The ducted propfan utilizes the blade design concepts of the unducted 
propfan, while incorporating a thin profile cowl which aids in maintaining the aero- 
dynamic loading in the blade tip region. An illustration of the aerodynamic charac- 
teristics associated with ducted propfans is given in Fig. 2.1. The (efficiency) and 
bypass ratio (BPR) of a ducted propfan are expected to lie somewhere in between 
the performance of current high bypass ratio turbofan engines (BPR 4-8) and recent 
unducted propfan demonstrator engines (BPR 25-40). In addition, ducted propfans 
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are expected to improve static thrust and acoustic signatures compared to unducted 
propfans, while permitting more flexible installation configurations. 

The advantages of the ducted propfan concept are offset by the additional drag 
imposed by the relatively large diameter cowl. Accurate prediction of the cowl drag 
has been the primary factor motivating the study of this concept in the CFD com- 
munity. A three-dimensional analysis for steady ducted propfan flows was recently 
presented by Hall and Delaney [9], and Hall et al. [10]. Williams et al. [11] utilized 
a frequency-domain panel method for predicting both steady and unsteady ducted 
propfan flows. Unfortunately, their analysis was limited to zero thickness airfoils and 
ducts, and cannot accurately account for the transonic, vortical nature of high speed 
ducted propfan flowfields. 

With the acceptance of CFD as a design analysis tool, the trend in CFD research 
has been to examine more demanding flow conditions and provide some insight into 
the complex flow physics which are often assumed away in the normal design analysis 
process. In the investigation of propfan aerodynamics, one area of advanced CFD 
research is aimed at gaining an understanding of the time-dependent velocity and 
pressure fields generated when the axis of the propeller is inclined at an angle of 
attack to the incoming freestream. While the freestream flow may be uniform, the flow 
relative to the rotating blades varies with circumferential position, and the resulting 
flow is inherently unsteady due to this nonaxial inflow. This is clearly an area of 
importance since aircraft engine propeller installations rarely operate under ideal 
(zero incidence angle) conditions. The periodic unsteady flow generated by the angle 
of attack operation may generate undesirable structural/aerodynamic interactions 
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and/or unacceptable acoustic noise levels, and therefore must be considered in the 
design process. Numerical methods for predicting the flow about unducted propfans 
at angle of attack have been demonstrated by Whitfield et al. [8] and Nallasamy and 
Groeneweg [12]. 

This document contains the Final Report for the ADPAC (Advanced Ducted 
Propfan Analysis Codes) 3D Euler/Navier-Stokes aerodynamic and grid generation 
analyses developed by the Allison Gas Turbine Division of the General Motors Corpo- 
ration under Task II of NASA Contract NAS3-25270. The objective of this study was 
to develop a three-dimensional time-dependent Euler/Navier-Stokes analysis for high- 
speed ducted propfan aircraft propulsion systems operating at angle of attack. This 
analysis consists of a grid generation scheme coupled with an advanced aerodynamic 
analysis code. The grid generation scheme is referred to as ADPAC-CHGRIDV2 or 
simply CHGRIDV2. The aerodynamic analysis is referred to as ADPAC-AOA or 
simply AOA. AOA utilizes a finite volume multiple-block four-stage Runge-Kutta 
numerical algorithm to predict the aerodynamics of ducted fan flows. Of particular 
interest was the ability to accurately predict the unsteady aerodynamics due to the 
angle of attack, and the complicated viscous flow interactions occurring between the 
rotating fan and the cowl. The use of a multiple grid block arrangement simplifies the 
calculation of the full rotor geometry required for angle of attack flows, and permits 
some unique grid arrangements for complicated ducted propfan geometries. 

Unducted propfans are analyzed using a single sheared H-type grid for each blade 
passage. The analysis for ducted propfans is based on a numerically coupled multiple- 
block grid arrangement with a body-centered C-type grid about the cowl, surrounded 



UNSTEADINESS RESULTING 
FROM INTERACTION WITH 
WINGS, STRUTS, OR 
COUNTER-ROTATION 


COWL BOUNDARY LAYER 


COWL WAKE 




SECONDARY. 
FLOW 



SPINNER/HUB/NACELLE 
BOUNDARY LAYER 


CLEARANCE 
TIP VORTEX 


UNSTEADINESS 
RESULTING FROM 
ANGLE OF ATTACK 



NONUNIFORM INFLOW 
RESULTING FROM 
INDUCED VELOCITIES 


BLADE SURFACE 
BOUNDARY LAYER 


BLADE WAKE 


Figure 2.1: Ducted propfan aerodynamic characteristics 
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by four H-type grid blocks for each blade passage. Predicted results using both grid 
systems were compared with available experimental data for several cases including 
a high-speed, high bypass 1.15 pressure ratio fan. 

To predict the flow about a ducted propfan at angle of attack using the analyses 
described in this document, the necessary sequence is: 

1. Define the geometry 

2. Generate a numerical grid for the domain of interest using CHGRIDV2. 

3. Run the Euler/Navier-Stokes code ADPAC-AOA to predict the steady aerody- 
namics. 

4. Rotate and duplicate the single-passage grid and steady state results using 
ADPAC-ROTCGRID and ADPAC-ROTCFLOW to provide full rotor initial 
data for the unsteady solution. 

5. Run the Euler/Navier-Stokes code ADPAC-AOA to predict the unsteady aero- 
dynamics. 

6. Plot and process the results as needed using ADPAC-AOAPLOT or other codes. 

Separate sections are provided in the chapters which follow to describe the basis 
and operation of the codes used in the steps above. A theoretical development of 
the grid generation algorithms for both ducted and unducted architectures is given 
in Chapter 2. The 3-D time-marching Euler/Navier-Stokes analysis is detailed in 
Chapter 3. In Chapter 4, a summary of the predicted results and verification studies 
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performed to validate the accuracy of the analysis is presented. A summary of the 
conclusions of this study is given in Chapter 5. 

It is worthwhile mentioning that the development and application of the codes 
described in this manual were performed on UNIX-based computers. All files are 
stored in machine-independent format. Small files utilize standard ASCII format, 
while larger files, which benefit from some type of binary storage, are written in 
a machine-independent format through the Scientific DataBase Library (SDBLIB) 
routines [13] The SDBLIB format utilizes machine- dependent input /output routines 
which permit machine independence of the binary data file. The SDBLIB routines 
are under development at the NASA Lewis Research Center. 

Most of the plotting and graphical postprocessing of the solutions was performed 
on graphics workstations. Presently, the PLOTSD [14], SURF [15], and FAST [16] 
graphics software packages developed at the NASA Ames Research Center are being 
extensively used for this purpose, and plot output has been tailored for this software. 
In addition, due to the increasing popularity of the PostScript page description lan- 
guage, and the variety of devices which can display PostScript - based output, a number 
of plotting procedures included in the ADPAC package utilize standard PostScript 


routines. 



9 


3. GRID GENERATION ALGORITHM 

In this chapter, the numerical algorithm forming the basis of the ducted/unducted 
propfan analysis grid generation scheme is described. The geometry and computa- 
tional domain are briefly described in the first section below. 

3.1 Computational Domain 

The problem of interest is either a ducted or unducted propfan geometry oper- 
ating at angle of attack, as shown in Fig 3.1. Geometric parameters are expressed in 
a cylindrical coordinate system referenced to the propfan axis. The axial coordinate 
lies along the propeller axis of rotation, the radial coordinate is perpendicular to the 
axis of rotation, and the circumferential coordinate sweeps in the counter-clockwise 
direction when viewed down the axis of rotation (i.e. looking downstream). This 
coordinate system is illustrated in Fig. 3.1. The hub and duct contours are assumed 
to be axisymmetric surfaces (no circumferential variation). This is presently a limita- 
tion of the grid generation scheme only. The aerodynamic analysis can be applied to 
nonaxisymmetric duct geometries presuming a suitable grid system can be generated. 
It is further assumed that the airfoil tips are fully enclosed within the duct for ducted 
configurations (although this limitation could be easily relaxed). For steady flow and 
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a periodic geometry, the computational domain may be reduced to a single blade 
passage, circumferentially, as shown in Fig. 3.1. The problem is further simplified by 
fixing the computational domain to the rotating blade elements. The resulting flow 
predictions are therefore based on the steady flow relative to the rotating blade. At 
angle of attack, however, we no longer expect the flow to be spatially perdiodic with 
respect to the blade pitch (circumferential spacing), and therefore the computational 
domain must be expanded to include the entire propfan. 

In order to implement a finite volume numerical solution of the unsteady aero- 
dynamics for a complete rotor system, the flowpath must be subdivided into a finite 
number of smaller elements within the overall region of interest. To enhance the qual- 
ity of the numerical solution, we further specify that these subdivisions be constructed 
m a relatively smooth and orderly manner with some arbitrary limitation on both 
the relative spacing between cells and the total number of cells in the computational 
domain. Each of these subdivisions will be used to define an elemental computational 
cell which will form the geometric basis of the finite volume aerodynamic solver to be 
described in a later chapter. 

The methods by which the computational cells are constructed differs slightly for 
ducted and unducted propfans. Unducted propfans utilize a single H-type grid per 
blade passage. This grid is then duplicated circumferentially to construct a multiple- 
grid block system for the complete rotor. Ducted propfans utilize a coupled system 
of five grid blocks per blade passage. An embedded C-type grid is wrapped about the 
leading edge of the cowl to enhance the resolution and accuracy of the aerodynamic 
predictions in this critical area. The C-type grid is surrounded by a network of 
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Ducted Propfan Analysis Computational Domain 



Figure 3.1: Ducted propfan analysis computational domain 
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four H-type grids which connect the inner C-grid to the outer computational domain 
boundaries. Separate sections are provided below to illustrate the details of the grid 
generation sequence for both unducted and ducted propfan geometries. 

3.2 Unducted Propfan Grid Generation 

The grid generation scheme for unducted propfans is based on a procedure 
originally developed by Mulac [17]. The approach here is to first construct a two- 
dimensional axisymmetric mesh in the meridional (z,r) plane, followed by a separate 
construction in the circumferential direction to determine the final three-dimensional 
mesh. 

The two-dimensional axisymmetric mesh is constructed by dividing the merid- 
ional plane into several subregions as shown in Fig. 3.2 for an unducted propfan. 
The regions of interest include the inlet, blade, exit, and outer flow regions. These 
regions are individually gridded in the two-dimensioned plane to satisfy the conditions 
of common grid points along region interface boundaries, and a suitable distribution 
of points within each region. 

The various regions differed in their individual construction. The distribution of 
points in any region is given axially and radially by one of three point distribution 
routines. These routines are described in the paragraphs below. All interpolations of 
coordinates were performed using spline interpolation routines. 

The distribution of the points on the meridional blade plane in both the radial 
and axial directions is determined by a simple geometric progression radiating from 
a symmetric centerline referred to as packing algorithm #1 described below. 



Unducted Propfan H-Grid Generation Subregions 



Figure 3.2: Meridional plane projection grid generation subregions for an unducted 
propeller 
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Packing algorithm #1 

This algorithm packs the mesh at both ends of the point distribution based only 
RAT. 

Given: 

D total projected length of cubic spline curve 
RAT ratio of adjacent cell projected lengths 

M number of points distributed across D (must be odd) 
the initial cell length of the progression is calculated as 


on 


DELT = 


D 

JL 


E|^/2) + l(JL4T)* 


The point distribution is then given by 


P(>) = E (DELTiRATii-V-, i = 2, (M/2) + 1 
1 = 2 


(3.1) 


(3.2) 


The remaining points (i = (M/2) + 2, M ) are determined from the known sym- 
metry of the point distribution. 


Packing algorithm =£2 


The second construction utilizes a slightly different packing algorithm. This 
algorithm packs the mesh either at the initial end of the point distribution (RAT > 
1) or at the terminal end (iL4T < 1). In this case, D, RAT , and DELT are specified 
along with N, the total number of curves in the region upon which the points are to 
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be distributed, leaving M to be determined. The value of M is calculated iteratively 
until the the following conditions are satisfied: 

/ n \ M 

{^ELf) rnax^ P} RAT)t 

2 = 1 
2=1 

where 

(' DELr)max = ma * (dELt) j=l,N 

(' DELf)min = ™ U ( DELT )j=l t N (3,3) 

These conditions guarantee that neither the adjacent cell ratios, nor their in- 
verses, exceed RAT. The resulting point distribution is then given by 

Pi}) = E ( DELT)(RAT)i , 2 = 1 , AT (3.4) 

i=2 

This approach is used axially in the inlet and exit regions, and radially in the 
outer boundary region. 


Packine algorithm #3 

The third packing algorithm is used to determine the axial point distribution 
between blade rows. This algorithm packs the mesh at both ends of the point distri- 
bution based on RAT and a given initial spacing at each end. The packing is denser 
at the ends for RAT > 1, and denser in the middle for RAT < 1. 

Given: 

D total projected length of cubic spline curve 
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RAT maximum ratio of adjacent cell projected lengths 
DELTTE initial cell width (forward blade row trailing edge) 

DELT LE initial cell width (aft blade row leading edge) 

N total number of curves upon which the points are to be distributed 
the following values are then calculated: 


DTE = 


DIE = 


DELTTE 

’deltteTdeltle 

DELTLE 

DELTTE + DELTLE 


D 

D 


( 3 . 5 ) 

( 3 . 6 ) 


The value of M is determined iteratively until the following conditions are satisfied: 


/ D \ M 

\DELt) max ~ S RAT% (3 ' 7) 

Z=I 

/ D \ M 

(: DELfLin £ S RAT ~' < 3 ' 8 ) 

Z = 1 

where: 

( D \ _ ( DTE \ ( DLE \ . 

\DELTJmax \DELTTEJ \DELTLE ) “ 1,N ^ 3 ’ 9 ^ 

(__£_\ _ . ( DTE \ ( DLE \ . 

\DELT) m in m * n \DELTTEJ \DELTLEJ ' 3 ~ lyN ^ 3 ‘ 10 ^ 

Again, this ensures that neither the adjacent cell ratios, nor their inverses exceed 

RAT . The point distribution is then given from the two ends of the curve by: 

P(i) = £ ( DELTTE){RAT)i , » = 1, M (3.11) 

3 = 1 


l 

P(2*M-l-l) = D- x: (DELTLE)(RATy, l = 1,M (3.12) 

j = 1 
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The blade region is generated initially. This region is defined by the axisymmetric 
projection of the blade hub and tip boundaries, and the blade leading and trailing 
edges. The number of points defining the blade region axially and radially is directly 
specified by the user. The point distributions are determined by packing algorithm 
#1. Once the blade points are determined, the inlet region pictured in Fig. 3.2 is 
computed next using packing algorithm #2. The exit region is constructed in the 
same manner as the inlet region. Finally, the radial distribution of points in the outer 
region is determined from packing algorithm #2. Once the meridional plane grid 
has been calculated, the full three-dimensional grid is constructed. The blade twist 
and thickness distributions are used to determine the circumferential coordinates for 
the blade regions. The remaining circumferential coordinates are based on utilizing 
packing algorithm #1 in the blade-to-blade plane. The circumferential variation of 
grid points upstream and downstream of the blade regions is adjusted to provide a 
smooth transition such that the grid lines become parallel to the axis of rotation away 
from the blade region. 

A sample grid for an unducted propfan based on this grid generation technique 
is given in Fig. 3.3. Once the grid for a single blade passage has been constructed, it 
is a simple matter to duplicate and rotate this grid circumferentially to complete the 
construction of the full rotor. The total number of gnd blocks in this case is IV, where 
N represents the number of rotor blades. Again, this clearly implies an axisymmetric 
hub surface, although the aerodynamic analysis does not require this limitation. 

The emphasis in this development, as in Mulac [17], was to maintain a reasonable 
grid quality, rather, than a specific number of grid points. Due to the algebraic 
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nature of the scheme, the construction is performed relatively easily and inexpensively. 
Unfortunately, this approach can have drawbacks under certain conditions. Since 
a portion of the grid generation sequence involves finding the intersection of two 
aribtrary curves, it has been observed that this can produce unexpected results when 
the two curves in question have nearly identical slopes near their intersection point. 
In some cases, it is possible to generate curves which have no intersection, and the 
grid generation scheme will clearly fail. In addition, blades or cowls with thick or 
blunt leading edges may not be well represented by this scheme when a relatively low 
number of grid points is used. Finally, the construction may occasionally have limits 
based on the number of grid points in a particular region due to the inability to satisfy 
one or more of the adjacent cell ratio constraints in the packing algorithms. In spite 
of these drawbacks, for most problems involving relatively conventional geometries, 
this approach has been found to be more than adequate. 

3.3 Ducted Propfan Grid Generation 

The grid generation scheme for ducted propfans utilizes a number of concepts 
developed in the previous section relating the construction of grids for unducted prop- 
fans. Again, the mesh for a single blade passage is generated by first constructing a 
two-dimensional grid in the meridional plane. For ducted propfans, a slightly differ- 
ent subdivision approach is utilized. Fig. 3.4 illustrates the boundaries of the various 
subregions used in the construction of the grids for ducted propfan geometries. The 
overall grid system consists of five separate grid blocks as shown by the thick lines 
defining the block boundaries. A single body-centered C-type grid is wrapped about 
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Figure 3.3: 


Sample grid for an unducted propfan 
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the leading edge of the cowl, represented by Block #3 in Fig. 3.4. The use of the 
body-centered grid is advantageous in that the grid skewness normally associated 
with simple sheared H-type grid systems n the vicinity of the cowl leading edge can 
be eliminated. The remaining four blocks (#1, #2, #4, and #5) are H-type grids 
which surround the C-type grid and connect the inner grid to the outer boundaries of 
the computational domain. Several of the subdomains of the individual grid blocks 
are similar to those used in the unducted propfan grid construction, so the extension 
of the algorithm for ducted propfans is obvious. The blade, mid-inlet, inlet, mid-exit, 
exit, and outer regions are all constructed in basically the same manner for the C-grid 
as they were for the H-grid. In this case, the axial locations of the upstream boundary 
of the mid-inlet and the downstream boundary of the mid-exit regions are specified 
as a fraction of the cowl axial chord upstream and downstream of the cowl leading 
and trailing edges, respectively. The outer boundary for the C-grid (block #3) and 
the upstream and downstream H-grids (blocks #2 and #4) is created from a line 
displaced from the outer surface of the cowl by an amount equal to the average of 
the blade leading edge and trailing edge tip gaps. This line is extended upstream and 
downstream at the radius of the ends. The near-cowl region is then discretized by the 
C-grid as shown in Fig. 3.4. The construction of blocks 1,2,4, and 5 follows directly 
from the development of the the subregion grids for unducted propfans, so no further 
discussion is presented here. Instead, the underlying algorithm for the C-grid (block 
3) is presented in the paragraphs below. 

Grid points for the C-grid block surrounding the cowl are determined in a three- 
step procedure. In the first step, the inner and outer boundary points are specified. 
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Ducted Propfan Multiple-Block C-Grid Subregions 



Figure 3.4. Axisymmetric plane projection gnd generation subregions for a ducted 
propeller 
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In the second step, the interior points are determined. In the third step, the interior 
grid points are clustered near the surface of the cowl to adequately resolve the cowl 
boundary layer for viscous flow calculations. Each step has several options each of 
which is described below. 

The outer boundary points of the C-grid are determined by the surrounding four 
blocks. This is based on the requirement of coincident points along all block bound- 
aries. The number of points along the upstream plane of the C-grid is determined by 
twice the number of radial grid lines in the C-grid (NPBCAB) plus one. The number 
of points from the cowl surface to the outer boundary of the C-grid (NPBCAB) is 
specified by the user, which also fixes the number of points along the C-grid block exit 
plane. The cowl surface points are determined through one of two procedures. The 
first method performs a simple interpolation of the input coordinates based on arc 
length around the cowl such that the relative arc length distributions of the cowl sur- 
face points and the corresponding outer boundary points are the same. Unfortunately, 
this can often lead to highly skewed mesh lines near the cowl surface. To circumvent 
this problem, an option was provided which allows the cowl surface grid points to 
“float” along the contour of the cowl in a manner that would impose orthogonality 
at the cowl surface to a line extending from the cowl point to the corresponding 
outer boundary point. The movement of the cowl surface grid points is controlled 
by a secant iteration procedure which optimizes the orthogonality of the surface grid 
point location as a function of arc length along the cowl, while maintaining a smooth 
transition through neighboring surface grid points. The new cowl coordinates are 
determined from the updated value of arc length through a linear interpolation of the 
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arc length and cowl coordinates originally specified. The resulting grid thus possesses 
improved orthogonality characteristics {dong the cowl contour, which is desirable in 
terms of solution accuracy. 

The secant iteration procedure described above for the cowl surface grid points 
is expressed as: 


1 _ s k , ( 0-0 

5 - =Si+ ""(^oM) 


( 3 . 13 ) 


where Sj is arc length measured clockwise around the cowl from a fixed reference 
location (1) to the point (i), k is the secant iteration count, and Oj is the measure of 
nonorthogonality. The orthogonality measure was based on the change of arc length 
Sgrid between the outer point and the corresponding inner point. The optimized 
point occurs when O j = dSg r ^d/dS ^ = 0. 

The starting values for the secant iteration are determined by using an initial 
point distribution obtained from the cowl surface grid point distribution scheme based 
on arc length (the first method described above). 

In order to avoid overlapping grid lines and to maintain stability, the new surface 
grid point locations were never allowed to migrate more than one third of the distance 
between the original point and the neighboring grid points. 


The interior grid points of the C-grid are determined through one of two meth- 
ods. The first method is a simple algebraic interpolation of coordinates between the 
inner and outer C-grid boundary point distributions. In the second method, the 
interior points are generated through the numerical solution of a set of elliptic equa- 
tions controlling a weighted distribution of grid smoothness, orthogonality, and grid 
point density based on a variational formulation originally developed by Brackbill 
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and Saltzman [18]. A brief description of this scheme is given below. 

For a two-dimensional grid, the following integral expressions may be derived to 
evaluate critical aspects of the overall grid quality in physical space: 


Grid smoothness: 

I* = J j l(V() 2 + (Vtj)^Jdzdr 

(3.14) 

Grid orthogonality: 

: ° = ff K V «) ' (Vn))izir 

(3.15) 

Grid point density: 

Iw = J J w(z,r)Jdzdr 

(3.16) 

where: 

N 

II 

P* 

II 

(3.17) 


J W't,) 

d(z,r) 

(3.18) 

and where V is the Cartesian gradient vector operator. The term w(z,r) 

is a user- 


specified function the magnitude of which is proportional to the desired grid point 
density in physical space. 

Obviously, the smoothest possible gnd is obtained when 1$ is minimized, the 
most orthogonal grid is obtained when I 0 is minimized, and the grid with the most 
desirable point density is obtained when Iyj is minimized. By minimizing a weighted 
sum of these terms, i.e: 

I = Zs ■+■ Colo + C W I W (3.19) 

the constants Co, C-w may be used to control the relative importance of orthogonality 
and point density, respectively, in the overall grid point distribution. By exchang- 
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ing dependent and independent variables, and applying the concepts of variational 
calculus for minimizing functions using the Euler-Lagrange equations, the following 
nonlinear coupled set of equations results: 


b l z tf + b 2 z fr + h z im + °l r & + a 2 r £i] + a 3 r w = (3.20) 

o 1 dw 

a l z tf + a 2 z (ij + a 3 z VV + c l r ££ + c 2 r £r) + c 3 r VV = ~ J (3.21) 

where the coefficients (o t *, 6j, Cj,t = 1, 3) are all functions of the coordinate derivatives 
as: 

a l = a sl + Coa-oi + Cwo- v i 
a 2 = a s2 + C<>a- 0 2 + Cwa v 2 
a 3 = a s3 + C<>a 0 3 + C w a v 3 
^1 = h 3l + Cob Q i + Cyjbyi 
b 2 = b s 2 + Cob 0 2 + Cwb v 2 
b 3 = b s 3 + Cob 0 3 + C w b v s 
C 1 = c sl + @oc 0 \ + Cwc-vl 
c 2 = c s 2 +'Coc 0 2 + Cwc-1,2 

c 3 = c s3 + ^o c o3 + Cw c v3 (3-22) 

where: 

a sl ~ — « a 2 = 2 ( aa )A <1,3 = —(00)7 

5,1 = (bb)a, b s 2 = —2(65)/?, 6,3 = (66)7 

c sl = C 00 ) 0 * c s2 = “ 2 (cc)0, C $3 = -(cc) 7 
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a ol = z y r yi 

a o2 = z ( r y + z y T £‘> 

a o3 = z t r £ 

b ol = z % 

b o2 = 2(2z^z^ + r^rrj), 

b oZ = z \ 

c ol = T 1 /» 

c o2 = 2 ( Z ( Z V + 2 r {r v ), 

c oZ = r \ 

a v\ = - z y r T]i 

a v2 = z ( r n + z V r (> 

a vz = ~ z a 

r 

II 

-3 tO 

b v2 = — 2r £ r y> 

Kz = r \ 

c vl — z yt 

c v2 = 2z £ s T)i 

c vZ = z \ 

(aa) = z£T£ + zrjrrj , (55) = + r 2 , 

(cc) = z 2 + z 2 

+ y%)/J z , P-- 

= (x£x v + y^/J 2 , 

7 = (x 2 + y 2 )/J' 


(3.23) 


When C 0 — 0 and C w = 0 the standard Laplace grid generation scheme results. 
This system is more complex than the usual Laplace-based grid generation schemes 
[19], but is still solvable using standard relaxation techniques. In this case, an iterative 
successive overrelaxation Gauss-Seidel solution technique is applied to solve the finite- 
difference equations resulting from a second-order central-difference approximation of 
the resulting equations. For example: 




*»+l, j ~ 2 Zjj + Xj_ hj 

(AO 2 


X (V 


x *+ hi+1 Z + g »-l J-l - **-1,7+1 

(A(A V ) 


(3.24) 

(3.25) 


Normally, the recommended approach is to set C 0 ,C W = 0 which effectively 
reduces the scheme to a Laplace solver. The additional complexity of the variational 
approach is presented here for completeness, but has not been found to be particularly 
useful for the present applications. 
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Once the interior points of the C-grid have been generated through the sequence 
just described, a reinterpolation of grid coordinates is performed along grid lines 
normal to the cowl surface to provide mesh clustering about the cowl for viscous 
flow calculations. The clustering is based on a one-sided Roberts transformation 
(see e.g. [20]). The initial grid distribution is used as the basis for the interpolation 
performed during the mesh clustering process. All interpolations are based on arc 
length S along a given grid line emanating from the cowl surface. Given the distri- 
bution of coordinates (z,r) along a given grid line, and the corresponding arc length 
distribution S, the clustered dsitribution of coordinates is determined as: 


c. = c. (/?+!)-(£ -!)[(£ +!)/(/? -PI 1 -; 

; ]max [08 + l)/(/J-l)]W + l 


(3.26) 


where is a user-specified parameter which controls the amount of clustering 
more clustering as j3 approaches 1), j = j / jmax where j is the numerical index 
of a particular point, and jmax is the maximum numerical index of the grid normal 
to the cowl surface. A corresponding reinterpolation of points is also performed for 
grid block #4 to maintain the continuity of grid points across the block boundary 
separating blocks #3 and #4. 

A sample multiple-block C-grid mesh network is shown in Fig. 3.5 for a ducted 
propfan geometry. 
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Figure 3.5: Sample multiple- block C-gnd for a ducted propfan 
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4. 3D EULER/NAVIER-STOKES NUMERICAL ALGORITHM 

This chapter contains a description of the time-dependent multiple-grid block 
Euler/Navier-Stokes ducted propfan aerodynamic analysis. The definitions of the 
pertinent variables used in this chapter may be found in the Nomenclature. 

4.1 Nondimensionalization 

To simplify the numerical treatment, all variables in the numerical solution are 
nondimensionalized by reference values as follows: 


z f 

z = T » r = r ’ 

L ref Kef 

* ^ 
«?* 
II 

ii 

£ 

v e = ** 

v ref 


P A* 

V — i P — > 

Pref Pref 

Cp 

' P= *r'f' 

H 

<3* 

Kef 


II 

2 s. 

*-s 

II 

* L ref 

w = — 

v ref 


(4.1) 


The reference quantities are defined as follows: 

L re j is the maximum diameter of the propfan blade 

Vref * s the frccstream relative total pressure 

P re f is the freestream relative total density 

a re f is determined from the freestream relative total conditions 
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— \JlPref / ~Pref 

v Te j is determined from the freestream acoustic velocity as 



/z re y is the freestream viscosity 
k re j is the freestream thermal conductivity 
R re f is the freestream gas constant 
T re j is the freestream temperature 

4.2 Governing Equations 

The numerical solution procedure is based on the strong conservation law form of 
the Navier-Stokes equations expressed in a cylindrical coordinate system. The Euler 
equations may be derived as a subset of the Navier-Stokes equations by neglecting 
viscous dissipation and thermal conductivity terms (i.e. - y and k = 0). By integrating 
the differential form of the Navier-Stokes equations over a rotating finite control 
volume, the following equations are obtained: 

/ U Q)iV + Li "<’ iQ) = / KdV + l viAQ) (4.2) 

where: 

L inv(Q) = J dA \hnv dA z + &inv dA r + (S inv - ruQ)dAg\ (4.3) 

and: 

LvisiQ) = + G v i s dA r + ByigdAg] 


(4.4) 
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The terms L^ nv and L v ^ a represent the cell face mass, momentum, and energy flux 
evaluations for the inviscid, and viscous components, respectively. 

The vector of dependent variables Q is defined as: 


Q = 


p 

pv z 

pv r 

P V B 

PH 


(4.5) 


where the velocity components v z ,v r , and vq are the absolute velocity components 
in the axial, radial, and circumferential directions relative to the propfan coordinate 
system, respectively (see e.g. - Fig. 3.1). The total energy function, e t , is defined 
as: 

* = (y-l)p + + V t + V b ( 4 -6) 

The individual flux functions are defined as: 


pvz 


pv r 


P v 0 

pVz+P 


pv z vr 


pv z vp 

pv z v r 

» ^inv — 

pv$+p 

> -®inv ~ 

pv T VQ 

rpv z VQ 


rpvrvg 


r{pv 0 + V ) 

. pv z H . 


. pv r H . 


. PVqH . 
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» ^vis — 

Trr 

» -®t ns ~ 
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T rB 
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- Qz . 


. Qr . 
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F = F(Q), G = G(Q), H = H(Q) 

F v = F V (Q), G v = G V (Q ), H v = H V (Q) 


( 4 . 8 ) 


( 4 . 9 ) 


The flux variables F, G, and H are determined at each grid cell interface by deter- 
mining the average (Q) of the cell-centered dependent variables from the individual 
finite volumes adjoining the interface. 

Finally, the cylindrical coordinate system source term is: 


K = 


0 

0 

PVfi+P 

r T ee 

0 

0 


( 4 . 10 ) 


It should be noted that in the numerical algorithm, the radius used in the cylindrical 
source term K is carefully formulated to guarantee numerical conservation for the 
radial momentum equation. That is, for a uniform stagnant flow, the radius in the 
radial momentum equation is chosen such that both sides of the radial momentum 
equation are equal. This ensures that small geometric errors do not corrupt the 
conservative nature of the numerical scheme. The total enthalpy, H t is related to the 
total energy by: 


h = h + z 

p 


( 4 . 11 ) 
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The viscous stress terms may be expressed as: 


T zz = 2fi^j^j+X v V-V, 

(4.12) 

--'[(£)*(&)]• 

(4.13) 


(4.14) 

rrr = 2 fi + A v V • V , 

(4.15) 

«-»[(SM£) -(?)]• 

(4.16) 

II 

to 

+ 

1 

+ 

e 

<1 

(4.17) 

, ,6T 

q z = v Z Tzz + VrT Z r + V $ T Z $ + 

(4.18) 

dT 

q r = VzTrz 4- v r T r r + VQT r Q + k - — , 

or 

(4.19) 

<10 = v z T 0 z + VrT 0 r + v e T ee + 

(4.20) 


where fi is the first coefficient of viscosity, \ v is the second coefficient of viscosity, 
and: 


V • V = — + ?2L + l du 9 + 

dz dr r d9 r 


(4.21) 


The remaining viscous stress terms are defined through the identities: 


t T z = T Z r, 

(4.22) 

*.»* 

II 

S- 

(4.23) 

ll 

* 

(4.24) 
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This integral form of the governing equations is applied to a generalized finite volume 
in physical space as shown in Fig. 4.1. The cell surface areas dA z , dA r , and dAg are 
calculated using the cross product of the diagonals of a cell face, and the cell volume 
is determined by a procedure outlined by Hung and Kordulla [21] for generalized 
nonorthogonal cells. 

In order to conveniently determine the viscous stress terms and thermal conduc- 
tion terms across an arbitrary cell interface, a generalized coordinate transformation 
is applied to the viscous stress terms as follows: 

t = r, 6), T) = i}(z, r, 0), ( = (( z , r, 0) (4.25) 

The chain rule may then be used to expand the various derivatives in the viscous 
stresses as: 

8^8 dij 8 8( 8 

dx dx 8£ dx drj dx (4.26) 

dr 8rd( drd-q^ drdC (4.27) 

— = Oil. ?1*L HiJL 

d$ ~ 80 8( + 80 drj + 80 8C ( 4,2S ) 

The transformed derivatives may now be easily calculated by differencing the vari- 
ables in computational space, and utilizing the appropriate identities for the metric 
differences (see e.g. [20]). 


4.3 Runge-Kutta Time Integration 

The time-stepping scheme used to advance the discretized equations is a four- 
stage Runge-Kutta integration. The solution proceeds as: 

I?1 = Q n - a^iHQ") + D(<?»)], 
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Figure 4.1: Three-dimensional finite volume cell 
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where: 


and: 


Q 2 = Q n - °2A<[£(<?i) + D(0 n )], 


<?3 = Q n - «3 a 'WQ2) + *>«?")], 


Qt=Q n - <* 4 A<[I(<? 3 ) + D(Q n )), 


<? n+1 = <?4 

(4.29) 

1 1 1 


8’ “ 2 = 4’ “3 = 2’ “4 = 1 

(4.30) 


HQ) = L inv(Q ) ~ L vis(Q) (4.31) 

Linear stability analysis indicates that this scheme is stable for all time incre- 
ments 6t which satisfy the stability criteria CFL < 2y/2. The CFL number may be 
defined in a one-dimensional manner as: 


CFL = 


At 




~£r 


(4.32) 


This factor is calculated for each coordinate direction, and then geometrically av- 
eraged to obtain the maximum allowable time increment for a given computational 
cell. 


For steady flows, an acceleration technique known as local time stepping is used 
to enhance convergence to the steady-state solution. Local time stepping utilizes the 
maximum allowable time increment at each point during the course of the solution. 
While this destroys the physical nature of the transient solution, the steady-state 
solution is unaffected and can be obtained more efficiently. For unsteady flow cal- 
culations, of course, a uniform value of the time step At must be used at every grid 
point to maintain the time-accuracy of the solution. 
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4.4 Fluid Properties 


The working fluid is assumed to be air acting as a perfect gas, thus the ideal gas 
equation of state has been used. Fluid properties such as specific heats, specific heat 
ratio, and Prandtl number are assumed to be constant. The fluid viscosity is either 
specified as a constant, or derived from the Sutherland (see e.g. [20]) formula: 


H = C 1 


3 

gg 

t + c 2 


( 4 . 33 ) 


The so-called second coefficient of viscosity \ v is fixed according to: 



( 4 . 34 ) 


The thermal conductivity is determined from the viscosity and the definition of the 
Prandtl number as: 


* = 2S 


( 4 . 35 ) 


4.5 Turbulence Model 


As a result of computer limitations regarding storage and execution speed, the 
effects of turbulence are introduced through an appropriate turbulence model and 
solutions are performed on a numerical grid designed to capture the macroscopic 
(rather than the microscopic) behavior of the flow. A relatively standard version of 
the Baldwin-Lomax [22] turbulence model was adopted for this analysis. This model 
is computationally efficient, and has been successfully applied to a wide range of 
geometries and flow conditions. 
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The effects of turbulence are introduced into the numerical scheme by utiliz- 
ing the Boussinesq approximation (see e.g. [20]), resulting in an effective calculation 
viscosity defined as: 

P effective = P laminar + A * turbulent (4.36) 

The simulation is therefore performed using an effective viscosity which combines the 
effects of the physical (laminar) viscosity and the effects of turbulence through the 
turbulence model and the turbulent viscosity turbulent • 

The Baldwin-Lomax model specifies that the turbulent viscosity be based on an 
inner and outer layer of the boundary layer flow region as: 

{ (^turbulent) inner i V — Vcrossover 

( \ ( 4 - 37 ) 

\Pturbulent) outer* V ^ Vcrossover 

where y is the normal distance to the nearest wall, and Vcrossover is the smallest 
value of y at which values from the inner and outer models are equal. The inner and 
outer model turbulent viscosities are defined as: 

(Pturb^inner = M (4.38) 

(Pturb) outer = KG cpP F wake F kleby (4.39) 

Here, the term / is the Van Driest damping factor 

1 = ky( 1 - e (~y + / A+ )) (4.40) 

u> is the vorticity magnitude, F wa ^ e is defined as: 


MwaJfee — ymaxFmax 


(4.41) 
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where the quantities ymaxi Fmax &re determined from the function 

F{y) = y\u\[l- e(~y + l A +)\ ( 4 . 42 ) 

The term y+ is defined as 

y 

The quantity F MAX is the maximum value of F(y) that occurs in a profile, and 
VMAX is the value of y at which it occurs. The determination of F MAX and VMAX 
is perhaps the most difficult aspect of this model for three-dimensional flows. The 
profile of F(y) versus y can have several local maximums, and it is often difficult to 
establish which values should be used. In this case, Fj^AX * s taken as the maximum 
value of F(y) between a y + value of 350.0 and 1000.0. The function F kieb is the 
Klebanoff intermittency factor given by 

F kleb(v) = I 1 + 5 > 5 ( ^ /e6y ) 6 ]— 1 (4.44) 

Umax 

and the remainder of the terms are constants defined as: 

A + = 26, 

Cep = 1 . 6 , 

C kleb = °* 3 » 
k = 0.4, 

K = 0.0168 (4.45) 

In practice, the turbulent viscosity is limited such that it never exceeds 1000.0 times 
the laminar viscosity. 
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The turbulent flow thermal conductivity term is also treated as the combination 
of a laminar and turbulent quantity as: 


^ effective ~ ^ laminar + ^ turbulent 


(4.46) 


For turbulent flows, the turbulent thermal conductivity k turbulent * s determined from 
a turbulent Prandtl number P T turbulent suc ^ 


Pr 


turbulent ~ 


C V ^turbulent 

L 

K turbulent 


(4.47) 


The turbulent Prandtl number is normally chosen to have a value of 0.9. 

In order to properly utilize this turbulence model, a fairly large number of grid 
cells must be present in the boundary layer flow region, and, perhaps of greater im- 
portance, the spacing of the first grid cell off of a wall should be small enough to 
accurately account for the inner “law of the wall” turbulent boundary layer profile 
region. Unfortunately, this constraint is typically not satisfied due to grid-induced 
problems or excessive computational costs, especially for time-dependent flow calcu- 
lations. A convenient technique to suppress this problem is the use of wall functions 
to replace the inner turbulent model function, and solve for the flow on a somewhat 
coarser grid. This technique has not been tested for the current application, but 
would appear to be a reasonable area for future research. Practical applications of 
the Baldwin-Lomax model for three-dimensional viscous flow must be made with the 
limitations of the model in mind. The Baldwin-Lomax model was designed for the 
prediction of wall bounded turbulent shear layers, and is not likely to be well suited 
for flows with massive separations or large vortical structures. There are, unfortu- 
nately, a number of applications for ducted and unducted propfans where this model 
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is likely to be invalid. This is also likely to be an area requiring improvement in the 
future. 


4.6 Artificial Dissipation 

An artificial dissipation operator ( D ) is added to the numerical scheme to control 
oscillations in the solution which result from the centered-difference approach in the 
flux derivative formulation. This problem is especially prevalent near shock waves, 
and it has been observed that the formulation of the dissipative term can have a 
significant influence on the final numerical solution. Jameson [23] demonstrated that 
a dissipative system combining second- and fourth- difference smoothing terms can ef- 
fectively eliminate undesirable numerical oscillations without destroying the accuracy 
of the solution. The scheme presented below is stable for all time steps satisfying the 
C Fi-related time step limitation 

CFL < 2\/2 (4.48) 

The dissipation operator is constructed in the following manner: 


D z (Q) = d i . , — d. i ., 

*+jv7,fc *— 


(4.49) 


d. x 

t+y,3,k 


where: 


{Ati) i+y,k 


{€2) i+fa,k AzQ i+y,k {ei) i+y,k A * Q i+y,k\ 


(4.50) 


K2max ^i+ij,bnj,k) 


(4.51) 
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( e 4) i+ l^ = -x(0,« -e i+ , jk ) 

v . _ ~ 2p i,j,k + Pz-lJ,fel 

,-7 ’ + P*-1J,&I 

Typical values for the second and fourth difference damping constants are: 


( 4 . 52 ) 

( 4 . 53 ) 


K 


2 


1 

4 



( 4 . 54 ) 


The term Afj represents a one- dimensional equivalent of the maximim allowable 
time step in the given coordinate direction. The use of this factor introduces an 
eigenvalue scaling into the dissipation operator which minimizes the added dissipation 
in coordinate directions which do not limit the stability of the algorithm. 

The damping scheme described above may be applied directly for inviscid flow 
calculations, but must be modified slightly for viscous flow calculations. As the mag- 
nitude of the physical viscous dissipation grows, the artificial dissipation is no longer 
required, and, in fact, can prevent convergence to the desired solution due to the 
complicated interaction between the physical and numerical dissipation structures. 
In order to deal with this problem, a controlling term is added to the artifical dissipa- 
tion operator to smoothly eliminate the damping term near solid walls. This scheme 
does not reduce the artificial damping in the freest ream, and could interfere with 
larger freestream vortical structures which may not require added dissipation. This 
implementation of the damping scheme is another area requiring future study. The 
dissipation in viscous regions is controlled by a simple one-dimensional function: 


(t2) i+l j,k = 

= Tmax(0, (k 4 ) - (e 2 ) 


i+T£ t j,k 


) 


( 4 . 55 ) 


( 4 . 56 ) 
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where: 

'V* • / ^ Tflf i 

T = . . ) ( 4 -57) 

Here m is the index of the coordinate direction along which the nearest wall is found, 
and m Umit * s 411 approximate grid reference location for the boundary layer edge. 
This simple model was found to be satisfactory for most calculations, and avoids the 
necessity for finding the true boundary layer edge (if one can actually be defined for 
complicated three-dimensional flows). 

The complete dissipation operator Z>— ^ is constructed as the sum of the dissi- 
pation operators in each of the respective coordinate directions as: 


Di,j,k ~ ( D z)ij,k + (®r)i j,k + ( D >)xj,k 


(4.58) 


4.7 Implicit Residual Smoothing 

Implicit residual smoothing is a technique commonly used for accelerating the 
convergence of explicit time-marching schemes applied to steady flow calculations. 
Since an unsteady flow calculation for a given geometry and grid is likely to be com- 
putationally more expensive than a similar steady flow calculation, it would be ad- 
vantageous to utilize this acceleration technique for time-dependent flow calculations 
as well. In recent calculations for two dimensional unsteady flows, Jorgensen and 
Chima [24] demonstrated that a variant of the implicit residual smoothing technique 
could be incorporated into a time-accurate explicit method to permit the use of larger 
calculation time increments without adversely affecting the results of the unsteady 
calculation. The implementation of this residual smoothing scheme reduced the CPU 
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time for their calculation by a factor of five. This so-called time-accurate implicit 
residual smoothing operator was then also demonstrated by Rao and Delaney [25] 
for a similar two-dimensional unsteady calculation. Although this “time-accurate” 
implicit residual smoothing scheme is not developed theoretically to accurately pro- 
vide the unsteady solution, it can be demonstrated that errors introduced through 
this residual smoothing process are very local in nature, and are generally not greater 
than the discretization error. 

In this study, a variant of the time-accurate implicit residual smoothing operator 
described by Jorgensen and Chima [24] was extended to three spatial dimensions 
and implemented in the unsteady solution procedure. The standard implicit residual 
smoothing operator can be written as: 

(1 - e z 6 zz )( 1 - e r £rr)(l - ie 0 S 00) R iJ t k = R i,j,k (4.59) 

where Rijfi (the residual) is expressed as: 

R ij,k = *WQij,k)inv + MQij,k)vis + £(<?;, j,k)\ (4-60) 

Here the differencing operator 8 is discretized as: 

SzzQ iJ,k = Qi+lJ,k - 2 Qij,k + Qi-l,j,k 

SrT Qi,j,k = ®ij+i,k - 2 ®ij,k + 

s eeQi,j,k = Qij,k+i ~ 2Q iJ,k + Qij,k - 1 

The reduction is applied sequentially in each coordinate direction as: 

. R ij,k = ( 1 ~ € * S zz)-l R ij,k 
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R Vj,h ~ (! ~ e ‘ Sr ')-l R lj,k 

*r/,k - e - n-iik 

= R iJ,k (4.61) 

where each of the first three steps above require the inversion of a scalar tridiagonal 
matrix. The residual smoothing operator is applied at the first and third stage during 
the four-stage Runge-Kutta algorithm. The time-marching scheme then becomes: 

Ql = Q n -a l R(Q n ) 

Q 2 = Q n - a 2 R(Ql) 

<?3 = Q" - «Z R (Q 2 ) 

Qi = Q n - “ 41 i(Q 3 ) 

<? B+1 = Ql (4.62) 

For steady flow calculations, a constant value of — gj. = — 2 is typically 

used to provide accelerated convergence. It can be shown that the Runge-Kutta time 
stepping scheme described in this report becomes unconditionally stable for any time 
step At when e satisfies: 


1 -• L 9 


(4.63) 


where e,- now varies throughout the grid, and CFL^ k represents the local value 
of the CFL number based on the calculation time increment At: • l, and CFL* . , 

l yJi K 

represents the maximum stable value of the CFL number permitted by the unmod- 
ified scheme. It is obvious then that the residual smoothing operator need only be 
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applied in those regions where the local CFL number exceeds the stability-limited 
value, and therefore the local variation of e,- • l is determined as: 


e ij,k 


= max 




CFL 


CFL* . , 

i,j,k 


M ^)2 _ : 


(4.64) 


In this approach, the residual operator coefficient becomes zero at points where the 
local CFL number is less than that required by stability, and the influence of the 
smoothing is only locally applied to those regions exceeding the stability limit. Prac- 
tical experience involving unsteady flow calculations suggests that for a constant time 


increment, the majority of the flowfleld utilizes CFL numbers less than the stability- 
limited value to maintain a reasonable level of accuracy. Local smoothing is therefore 
typically required only in regions of small grid spacing, where the stability-limited 
time step is very small. Numerical tests both with and without the time-accurate 
implicit residual smoothing operator for the flows of interest in this study were found 
to produce essentially identical results, while the time-accurate residual smoothing 
resulted in a decrease in CPU time by a factor of 2-3. 


4.8 Boundary Conditions 


Inflow and exit boundary conditions are applied numerically using characteristic 
theory. A one- dimensional isentropic system of equations is utilized to derive the 
following characteristic equations at an axial inflow/outflow boundary: 


de- 

bt 


< \ dC ~ n 

- (Vz - a)— = 0, 


dC+ . .3C+ 

1- (vz + a)-T — 


= 0 


(4.65) 


dt 


dz 


(4.66) 
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where: 

2 a i 2a 

C = v* r, C+ = v z + (4.67) 

7-1 7-1 ' 

In order to efficiently process boundary information in the numerical solution, phan- 
tom cells are located just outside the computational domain to permit the unmodified 
application of the interior point scheme at near boundary cells. Boundary condition 
information is effectively introduced into the solution by properly controlling the 
dependent variables in the phantom cells while permitting the application of the 
standard interior point scheme at near boundary cells. 

For subsonic normal inflow, the upstream running invariant C~ is extrapolated 
to the inlet, and along with the equation of state, specified total pressure, total 
temperature, and flow angle (used to specify the Jingle of attack), the flow variables 
at the boundary may be determined. It should be mentioned that the effective inflow 
angle may vary for a given block as it rotates about the axis, and therefore the inflow 
angle is actually a function of circumferential position, 9. At the exit, a static pressure 
is specified at the hub for internal flows, and at the outer boundary for external flows. 
The remaining pressures along the outflow boundary are calculated by integrating the 
radial momentum equation: 

dp pv# 

8t r 

In this case, the downstream running invariant C+ is used to update the phantom cells 
at the exit boundary. Far-field boundaries also use this characteristic technique based 
on whether the local flow normal to the boundary passes into or out of the domain. 
The solid surfaces (hub, cowl, airfoils) must satisfy flow tangency for inviscid flow: 

F • n = 0 



( 4 . 69 ) 
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or no slip for viscous flows: 

v z = 0, v r = 0, VQ = ru> (4.70) 

In both cases, we specify no convective flux through the boundary (an impermeable 
surface), and hence, only pressure is needed at the phantom cell. The pressure may 
be extrapolated, or updated using a variant of the normal momentum equation. In 
this case, extrapolation was found to be the most effective technique based on rapid 
convergence and adequate results. In addition, solid surfaces are also assumed to be 
adiabatic, which implies that the normal temperature gradient is also zero. 

4.9 Multiple-Block Coupling 

For the multiple-block C-grid scheme, the solution is performed on a single grid 
block at a time. Special boundary conditions along block boundaries are therefore 
required to provide some transport of information between blocks. This transport 
is provided through a simple procedure which relies on the fact that the grid block 
boundaries have coincident grid points. Since the standard interior point scheme 
is performed at the first grid cell off the block boundary, all that is required is to 
determine the fluxes along the block boundaries themselves. Phantom points are 
provided along each side of a block boundary to accommodate the flux calculation at 
the block boundary cell face. The flow variables in the phantom cells are provided 
by interrogating the corresponding cell in the adjacent block when a particular value 
is required. After each stage of the Runge-Kutta integration for each block, the 
block boundary phantom cells are updated by utilizing the new dependent variable 
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values from the adjacent blocks. This direct specification technique was compared 
to a characteristic-based method similar to the inlet and exit boundary condition 
routines, and was found to provide enhanced convergence behavior. 

For ducted configurations, the grid blocks used to discretize the blades are rotat- 
ing, while the grid blocks discretizing the cowl and outer flow regions are stationary. 
Under these circumstances, flow variables for the phantom points are found through 
a circumferential interpolation of the neighboring computational blocks based on the 
known angular position of the blades. This technique presumes that the meridional 
coordinates of the grid points along the block boundaries are coincident. It should 
be mentioned that the interpolation technique does not strictly enforce a global con- 
servation across the block boundary, although this has not posed a problem in the 
calculations performed to date. 

Artificial damping is applied at the block boundaries by neglecting the fourth 
order derivative term due to the lack of a four point differencing stencil at the bound- 
aries. Implicit residual smoothing is applied at the block boundary by imposing a 
zero residual gradient (i.e. (dR/dz) = 0.0) condition at the boundary. 

4.10 Solution Procedure 

The numerical solution is performed in an identical manner for both unducted 
and ducted propfan geometries. Assuming that the numerical grid and flow parame- 
ters are known, the time-marching procedure may begin from some set of initial data. 
This initial data is specified as a uniform flow, or may be introduced from a previous 
solution. The time-marching procedure is applied iteratively to update the flow vari- 
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ables as the solution proceeds. Steady state solutions are deemed converged when 
the average residual R has been reduced by a factor of 10~ 3 , or when the residual 
has stopped converging. 

For unsteady flows, it is usually best to first obtain a steady state solution for a 
single blade passage, and then duplicate this flowfield for the remaining blade passages 
to construct the initial data for the full rotor for an unsteady solution. The unsteady 
solution may then be advanced in time with a uniformly specified time step, until a 

time-periodic solution is achieved. This normally requires two complete revolutions 
of the rotor. 



5. RESULTS 


Several numerical results from the multiple-block 3D Euler/Navier-Stokes anal- 
ysis code AO A described in Chapter 4 are presented in the sections which follow. In 
most cases, steady flow results are given to verify the accuracy of the formulation be- 
fore presenting the unsteady computations. Viscous and inviscid flow calculations are 
presented initially for a 2-bladed SR7 propfan for which both steady state and time- 
dependent experimental airfoil surface static pressure distribution data were available. 
Next, the Euler algorithm is utilized to predict the flow about the complete 8-bladed 
SR7 propfan geometry, illustrating the time-dependent formation and dissipation of 
a blade passage shock system. A demonstration of unsteady inviscid flow about 
a generic ducted propfan geometry at angle of attack is presented in the following 
section. Finally, results from viscous flow calculations are compared extensively with 
high speed experimental data for a 1.15 pressure ratio ducted fan configuration. Each 
case is discussed separately in the sections which follow. A summary of the compu- 
tational statistics (CPU time, grid size, etc.) associated with each of the calculations 
presented in this report is given in Table 5.1. 
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Table 5.1: Summary of computational characteristics for ADPAC test cases 


Test Case 

Steady/ 

Unsteady 

Viscous/ 

Inviscid 

# Blocks 

Grid Size 

Mach # 

Cray-II CPU 
(minutes) 

Modane 

Steady 

Inviscid 

1 

135,000 

0.2 

48 

Modane 

Steady 

Viscous 

1 

135,501 

0.2 

153 

Modane 

Unsteady 

Inviscid 

2 

15,180 

0.5 

54 

(Coarse) 

Modane 

Unsteady 

Inviscid 

2 

195,734 

0.5 

2322 

(Fine) 

Modane 

Unsteady 

Viscous 

2 

131,936 

0.5 

2546 

SR7 

Unsteady 

Inviscid 

8 

237,160 

0.8 

2136 

Ducted SR7 

Unsteady 

Inviscid 

40 

450,000 

0.8 

2640 

NASA Fan 

Steady 

Viscous 

5 

183,414 

0.75-0.85 

210 

NASA Fan 

Unsteady 

Viscous 

60 

488,320 

0.20 

3361 


5.1 SR7 2-Bladed Propfan Modane Tests 


In order to verify the numerics of the aerodynamic solver, several initial calcu- 
lations were performed for steady and unsteady flows about unducted propfans. A 
number of experimental studies have been performed for the unducted SR7 propfan 
geometry ranging from scaled wind tunnel tests to in-flight measurements. The SR7 
propfan design utilizes 8 blades with 41 degrees of sweep at the tip. A special hub 
contour is employed to eliminate flow choking at the hub. A description of the SR7 
design parameters is given in Fig. 5.1. 

The first extensive comparison of predicted results for steady and unsteady un- 
ducted propfan flows was performed based on the test of the SR7 airfoil at the Modane 
wind tunnel test facility reported by Bushnell et al. [26]-[27]. In the Modane tests, 
the propfan driver did not have enough power to drive the full 8-bladed configuration, 
so a 2-bladed version was tested instead. The airfoil surfaces were instrumented to 



Characteristic 


Number of blades 

lip sweep angle, a deg 

Mode' diameter, cm (In.) 

Tip speed, m/sec (ft/sec) (at 10.669 m(35 000 ft) 

Power loading, kW/m? (shp/ft?) (i.S.A. altitude 
Activity factor, AF ' 

! integrated design lift coef f Iclent ,c M 
Airfoils . NASA 

Ratio of nacelle maximum diameter to propeller diameter 
Cruise design Mach number 
Cruise design advance ratio 
Cruise design power coefficient 


62.23 (; 
243.0 ( 
256.85 


16 and 


a Geometr1c measurement frrom planform. 


Figure 5.1: SR7 propfan design characteristics 
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permit the measurement of steady and unsteady airfoil surface static pressures. 

Both steady and unsteady, viscous and inviscid flow calculations were performed 
for this geometry to facilitate a comparison with the experimental data. An illus- 
tration of the numerical grid used for the steady state inviscid and viscous flow cal- 
culations are presented in Figs. 5.2 and 5.3, repectively. The inviscid grid contained 
120x45x25 points in the axial, radial and circumferential directions, repsectively, while 
the viscous grid numbered 93x47x31. The viscous grid was slightly smaller in the axial 
direction to reduce overall CPU time. The clustering of grid lines near solid surfaces, 
necessary to resolve the blade surface boundary layer flow, is clearly evident in the 
viscous grid. 

Steady flow predictions were compared with experimental data for a nominal 
power coefficient of 0.250. The freestream Mach number was 0.2, with an advance 
ratio of 0.881. This particular flow condition gives rise to a strong leading edge vortex 
which dominates the flow characteristics of the suction surface of the airfoil, and was 
therefore considered a challenging test case for numerical simulation. 

Predicted steady flow (zero angle of attack) airfoil surface static pressure coeffi- 
cient distributions for both inviscid and viscous flow are compared with experimental 
data at spanwise locations of 28.4% and 93.1% ( r/R ) in Figs. 5.4-5.5, respectively. 
Both Euler and Navier-Stokes solutions produced a leading edge vortex, evidenced 
by the large pressure drop near the leading edge of the suction surface in Figure 5.4. 
The viscous flow calculation more accurately reproduces the experimental findings 
in this region, presumably due to a more accurate representation of the leading edge 
vortex, and in part due to the blockage imparted by the airfoil surface boundary layer. 
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Further comparisons with experimental data indicate that the exact location of the 
origin of the leading edge vortex is not adequately resolved in either calculation, and 
indicates that spanwise grid refinement can play a larger role in the accuracy of the 
results than either circumferential or axial grid refinement. The predicted trajectory 
of the leading edge vortex was also somewhat different than the trajectory indicated 
by the experimental data. This feature is best illustrated in the comparison of pre- 
dicted airfoil suction surface static to total pressure ratio contours given in Fig. 5.6. 
The predicted contours display a remarkable similarity to the experimental distribu- 
tion. The leading edge vortex path is indicated by the region of low pressure bending 
across the suction surface from approximately 30% span towards the airfoil tip. The 
experimental contours suggest that the leading edge vortex remains attached to the 
leading edge, while both viscous and inviscid contours suggest that the vortex drifts 
further aft axially on the blade. This slight discrepancy may be due to several factors. 
The exact deflected airfoil shape and blade setting angle are unknown (experimental 
uncertainty is on the order of 1 degree), which introduces a significant source of error 
in the calculation. In addition, to truly capture the detailed physics of this compli- 
cated vortical flow would require grid densities in excess of the present grid system, 
and therefore grid resolution must be considered a potential source of error. 

The viscous flow predictions for this case show some striking aerodynamic fea- 
tures resulting from the complicated flow over this highly three-dimensional airfoil. 
The leading edge vortex is clearly depicted in the particle flow trajectories illustrated 
in Fig. 5.7. The influence of the leading edge vortex is depicted by the color static 
pressure contours used to define the blade in this figure. During the development of 
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Figure 5.2: 2-bladed SR7 propfan geometry and steady state inviscid calculation 

grid for Modane test comparison 
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Figure 5.3: 2-bladed SR7 propfan geometry and steady state viscous calculation grid 
for Modane test comparison 
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Figure 5.4: Comparison of predicted and experimental airfoil surface static pressure 

coefficient distributions for 2-bladed SR7 propfan Modane test (28.4% 
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Figure 5.5: Comparison of predicted and experimental airfoil surface static pressure 

coefficient distributions for 2-bladed SR7 propfan Modane test (93.1% 
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Leading Edge Vortex 
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Figure 5.7: Illustration of predicted turbulent flow leading edge vortex particle tra- 

jectory traces for for 2-bladed SR7 propfan 
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the viscous flow solution scheme, both laminar and turbulent viscous flow predictions 
were performed for this case, thus allowing some insight into the effects of turbulence 
for highly loaded propfan airfoils. An interesting depiction of the results of this study 
is found in the suction surface shear flow patterns for both laminar and turbulent 
flow predictions given in Fig. 5.8. The shear flow patterns were generated by releas- 
ing particles at the first grid point off the surface and tracing the near surface flow 
pattern by restricting the particle paths to the first grid plane off the suction surface. 
This technique is roughly equivalent to a surface oil flow visualization technique. The 
suction surface flow patterns for both predictions clearly reveal the radial secondary 
flows generated by the leading edge and tip vortices. In addition, the laminar cal- 
culation reveals a small leading edge separation region near the propfan root leading 
edge. The turbulent flow results shows no such phenomena. As expected, the lam- 
inar result also reveals a large region of separated flow near the trailing edge of the 
blade. The low momentum fluid in the separation zone migrates radially outward, 
causing a complicated interaction of leading edge vortex, tip vortex, and radial migra- 
tion flows near the tip. By comparison, the turbulent flow prediction shows a much 
smaller region of separation at the trailing edge, and less radial migration, although 
the interacting flows near the tip are no less complicated. These results demonstrate 
the detailed flow patterns which may be captured with advanced, three-dimensional 
viscous flow calculations. 

All of the results for the SR7 presented thus far have been based on a single 
blade passage grid solution. Unsteady flow predictions for the SR7 2-bladed propfan 
were initiated from steady flow results similar to those already presented, although 
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slightly fewer grid points were utilized to minimize the computational cost. For the 
unsteady predictions, each single passage grid system was rotated and duplicated to 
discretize the adjacent blade passage, and the full features of the multiple-grid block 
AOA solver were utilized to track the solution for the complete rotor system. For the 
time-dependent calculations, the flow Mach number was 0.5, and the angle of attack 
was 3.0 degrees. Time-dependent inviscid flow predictions were generated for both 
a coarse and a fine grid to illustrate the effects of grid density on the unsteady flow 
predictions. The coarse grid utilized 46x15x11 grid points per blade passage, while 
the fine grid utilized 77x41x31 points per blade passage. Viscous flow predictions 
were performed for a single grid utilizing 56x38x31 points per blade passage. 

A steady state solution was used to initialize the time-dependent calculation. 
For reference, comparisons of predicted and experimental steady flow (zero angle of 
attack) airfoil surface static/inlet total pressure ratio distributions for this case are 
given at 28.4% and 94.4% span in Figs. 5.9-5.10, respectively. As in the previous 
steady flow comparisons, the agreement between predictions and experiment is good. 

Both viscous and inviscid solutions were then initated from the corresponding 
steady state (zero angle of attack) solutions. Roughly two complete revolutions of 
the rotor were required to reach a time-periodic solution. The results presented here 
correspond to the third revolution. 

Several cycles of both the coarse and fine grid inviscid as well as the fine grid 
viscous predicted phase-resolved time- dependent pressure histories are compared with 
experimental data at 64% radial span and 4.9% and 36.7% chord on the airfoil suction 
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Figure 5.9: Comparison of predicted and experimental airfoil surface static/total 

pressure ratio distributions for 2-bladed SR7 propfan (28.4% span, 
M=0.5, J=3.06). 
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Figure 5.10: Comparison of predicted and experimental airfoil surface static/total 

pressure ratio distributions for 2-bladed SR7 propfan (94.4% span 
M=0.5, J=3.06). ’ 
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(camber) surface in Figs. 5.11-5.12, respectively. The inviscid coarse grid prediction 
is observed to be reasonably accurate in the region downstream of the blade leading 
edge, but fails to correctly predict the amplitude of the unsteady pressure fluctua- 
tions in the vicinity of the leading edge. This is clearly attributable to insufficient 
grid resolution, as both the viscous and inviscid fine grid results more closely match 
the experimental data. This trend is evident in nearly every data comparison for this 
case. A similar set of plots is given at 64% span and 4.9%, 10.0%, and 63.3% chord 
for the airfoil pressure (face) surface in Figs. 5.13-5.15, respectively. The comparison 
of time-dependent pressure histories at 4.9% span in Fig. 5.13 is disturbing, since 
the fine grid solutions show a larger deviation from the experimental measurements 
than the coarse grid results. This deviation becomes even more perplexing due to the 
excellent agreement with experimental data at 10.0% chord given in Fig. 5.14. Based 
on these observations, it was concluded that the deviation in Fig. 5.13 was due either 
to experimental error, or is a result of some flow phenomena which is not adequately 
simulated (such as a flow separation or time-dependent vortex development) in the 
present solution strategy. Further comparisons of predicted and experimental data 
for the two-blade SR7 propfan suction surface at 91% span and 27.9% and 69.8% 
chord are given in Figs. 5.16-5.17, respectively. Again, the fine grid results show very 
good agreement with the experimental data, especially near the leading edge, where 
pressure fluctuations are at a maximum. Finally, a similar comparison is given for 
the pressure surface at 91% span at 27.9% and 89.8% chord in Figs. 5.18-5.19, re- 
spectively. It was observed that the predicted results demonstrate a slight phase shift 
compared to the experimental data which is accentuated in the comparisons at the 
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outer (91%) spanwise locations. It appears that the use of the fine grid in the nu- 
merical calculation reduces this phase shift in most cases. The difference between the 
experimental and calculation blade shapes under aerodynamic and rotational loading 
could also be a factor in this discrepency. An illustration of the three-dimensional 
viscous flow resulting for this case is given in Fig. 5.20. This figure illustrates the 
2-bladed propfan shaded by instantaneous static pressure color contours and particle 
trajectories emanating from the blade tips. The asymmetric aerodynamic loading 
induced on the propfan spinner and hub by the angle of attack is clearly evident, as 
is the upsweep of the instantaneous particle trajectories from the blade tips. The 
viscous flow results for this case did not differ significantly from the fine grid inviscid 
flow results, presumably because no significant viscous flow phenomena (massive flow 
separation, vortical development, etc.) axe present during the unsteady cycle at this 
operating condition. A more interesting comparison would likely result for larger an- 
gles of attack or increased blade loading levels, both of which will be areas for future 
investigation. 


5.2 SR7 8-Bladed Propfan 

A second unsteady calculation for the SR7 propfan blade geometry was performed 
for the full 8-bladed propfan configuration. Steady state inviscid calculations based 
on this geometry were performed over a wide range of advance ratios and blade setting 
angles for flight Mach numbers of 0.7 and 0.8 under Task I of this contract [10]. 

This test, based on the SR7 8-bladed propfan geometry, corresponds to the nu- 
merical test results presented by Nallasamy and Groeneweg [12]. Nallasamy and 


Static Pressure Difference (P-Psteady)/PO 



Figure 5.11: Comparison of predicted and experimental airfoil surface time- 

dependent static/total pressure ratio history for 2-bladed SR7 propfan 
(suction side, 4.9% chord, 64.0% span, 3 degrees angle of attack) 




Static Pressure Difference (P-Psteady)/PO 


0.150 -i 


0.100H 


Experiment 

ADPAC-Euler (Fine Grid) 

ADPAC-Euler (Coarse Grid) 

ADPAC-Navier-Stokes (Viscous Grid) 


0.050 — I 


0.000 H 


-0.050H 



- 0 . 100 - 


I I 1 

180.00 540.00 900.00 

Nondimensional Time or Rotation 


-180.00 


1 

1260.00 


1620.00 


Figure 5.12: Comparison of predicted and experimental airfoil surface time- 

dependent static/total pressure ratio history for 2-bladed SR7 propfan 
(suction side, 36.7% chord, 64.0% span, 3 degrees angle of attack) 



Static Pressure Difference (P-Psteady)/PO 


0.150 


Experiment 

ADPAC-Eiiler (Fine Grid) 

ADPAC-Euler (Coarse Grid) 

ADPAC-Navier-Stokes (Viscous Grid) 


0.050 


0.000 H 


- 0 . 050 H 


- 0.100 



- 180.00 


i T 

180.00 540.00 900.00 

Nondimensional Time or Rotation 


1260.00 


Figure 5.13: Comparison of predicted and experimental airfoil surface time- 

dependent static/total pressure ratio history for 2-bladed SR7 propfan 
(pressure side, 4.9% chord, 64.0% span, 3 degrees angle of attack) 
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Figure 5.14: Comparison of predicted and experimental airfoil surface time- 

dependent static/total pressure ratio history for 2-bladed SR7 propfan 
(pressure side, 10.0% chord, 64.0% span, 3 degrees angle of attack) 
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Figure 5.15: Comparison of predicted and experimental airfoil surface time- 

dependent static/total pressure ratio history for 2-bladed SR7 propfan 
(pressure side, 63.3% chord, 64.0% span, 3 degrees angle of attack) 
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Figure 5.16: Comparison of predicted and experimental airfoil surface time- 

dependent static/total pressure ratio history for 2-bladed SR7 propfan 
(suction side, 27.9% chord, 91.0% span, 3 degrees angle of attack) 
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Figure 5.17: Comparison of predicted and experimental airfoil surface time- 

dependent static/total pressure ratio history for 2-bladed SR7 propfan 
(suction side, 69.8% chord, 91.0% span, 3 degrees angle of attack) 
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Figure 5.18: Comparison of predicted and experimental airfoil surface time- 

dependent static/total pressure ratio history for 2-bladed SR7 propfan 
(pressure side, 27.9% chord, 91.0% span, 3 degrees angle of attack) 
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Figure 5.19: Comparison of predicted and experimental airfoil surface time- 

dependent static/total pressure ratio history for 2-bladed SR7 propfan 
(pressure side, 89.8% chord, 91.0% span, 3 degrees angle of attack) 
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Figure 5.20: Instantaneous propfan surface static/total pressure ratio contours and 

blade tip particle trajectories for 2-bladed SR7 propfan at angle of 
attack (M=0.5, J=3.06, angle of attack = 3 degrees). 
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Groeneweg predicted the unsteady inviscid flow about an 8-bladed SR7 propfan at a 
Mach number of 0.8, advance ratio of 3.12, and 4.6 degrees angle of attack. The 3/4 
radius blade setting angle was 60.2 degrees. Their predictions for this case indicated 
the time-dependent formation of a passage shock which develops during the down- 
ward (into the angled freestream) motion of the propeller. The shock wave ultimately 
moves forward through the blade passage, and then dissipates as the propeller begins 
to swing upward (away from the angled freestream). Although no experimental data 
were available for comparison, the unsteady shock motion and previous predictions 
make this an interesting case to study. 

An unsteady inviscid calculation was performed for this unducted propfan on a 
system of numerical grids containing 8 blocks (one for each blade passage). The grid 
system utilized a total of 237,160 points as shown in Fig. 5.21. Instantaneous static 
pressure contour plots of a blade-to-blade grid surface corresponding to roughly 50% 
span are given for each of the eight blade passages consecutively in the direction of 
rotation in Figs. 5.22(a)- 5.22(h), respectively for the inviscid flow prediction. These 
plots serve to illustrate the differences in blade passage pressure during each blade 
pitch rotation of the propfan. The passages begin about one blade pitch before top 
dead center and continue in the direction of rotation into the downward sweep. The 
convergence of the contours in Figs. 5.22(a)- 5.22(c) clearly illustrates the formation 
of the passage shock. Figs. 5.22(d)- 5.22(f) portray the motion and simultaneous 
dissipation of the shock as it moves forward through the passage, until it ultimately 
disappears and renews the time-dependent cycle (Fig. 5.22(h)). This sequence is 
qualitatively identical to the observations reported by Nallasamy and Groeneweg [12] 


e-a 
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in their calculations. An instantaneous depiction of the predicted propfan surface 
static/total pressure ratio contours is given in Fig 5.23. The asymmetric hub loading 
and blade passage shock movement are all easily identifiable in the color contours. Of 
particular interest is the noticable change in blade leading edge stagnation point at the 
airfoil root due to the angled freestream. It would appear that the root airfoil section 
should be particularly insensitive to incidence angle for good off-design performance. 

A rotational history of the predicted single blade power coefficient history is given 
in Fig. 5.24(a). A reproduction of a similar plot from Ref. [12] is also illustrated in 
Fig. 5.24(b) for comparison. The present prediction is in good qualitative agreement 
with Nallasamy and Groeneweg’s results. Quantitatively, the ADPAC results indicate 
a slightly larger peak to peak variation in the single blade power coefficient history. 
This discrepency is more than likely due to a combination of effects resulting from 
differencies in grid resolution, numerical damping, and/or time step size between the 
two calculations. 


5.3 Ducted Propfan Test Case 


Since the ducted propfan is a relatively new area of investigation, there are few 
available experimental or numerical data upon which the present algorithm may be 
compared. Calculations based on existing fan geometries do little to illustrate the 
aerodynamics of advanced ducted propfan geometries, and therefore an unsteady cal- 
culation based on a fictitious geometry is initially examined here. The test geometry 
was developed based on the SR7 blade geometry. Since propfan design philosophy 
utilizes cascade flow concepts for the root airfoil sections, it is unlikely that the ad- 



Figure 5.21: Full rotor grid for SR7 time-dependent inviscid calculations 
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Figure 5.22: 


Predicted instantaneous propfan inviscid blade passage static pressure 
contours for 8-bladed SR7 propfan at angle of attack (M=0.8) 
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Figure 5.23: Predicted instantaneous propfan surface static/total pressure ratio con- 

tours for 8-bladed SR7 propfan at angle of attack (M=0.8, angle of 
attack = 4.6 degrees) 
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Figure 5.24: Predicted rotational single blade power coefficient histories for 8-bladed 
SR7 propfan at angle of attack (M=0.8, angle of attack = 4.6 degrees) 
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dition of an outer duct would significantly effect the aerodynamics in this region, so 
this geometry was unaltered. However, the near tip region of the blade is not likely to 
be adequately suited for ducting, and, therefore, the blade geometry was created by 
discarding the outer 20% of the SR7 airfoil, and enclosing this clipped eight-bladed 
propfan in a low profile duct. The length of the duct was chosen to be representa- 
tive of a high bypass ratio unmixed exhaust flow ducted propfan configuration. This 
geometry is conceptually equivalent to the types of designs expected for the ducted 
propfan concept. 

Inviscid aerodynamic calculations were again performed for both a coarse grid 
and a fine grid solution for the ducted SR7 geometry at a Mach number of 0.8, advance 
ratio of 3.06, and 5 degrees angle of attack. Again, the 3/4 radius blade setting angle 
was 60.2 degrees. Only the results from the fine grid solution are presented here. 

The fine grid solution employed 450,000 grid points and is illustrated in Fig. 5.25. 
For this configuration, the mesh was comprised of 40 separate grid blocks (5 per 
blade passage). A series of instantaneous blade passage static pressure contour plots 
at roughly 50% span are given in Figs. 5.26(a)- 5.26(h) for the fine grid calculation. 
It is immediately obvious that a rather strong passage shock has formed as a result 
of these flow conditions. It is apparent that the blade setting angle is not optimally 
adjusted for this case. The shock oscillates slightly during the course of the revolution, 
but ultimately shows little change as a result of the rotation, and the flow remains 
choked in each passage throughout the cycle. This is a result, in part, of the flow 
straightening effect induced by the duct, which aids to minimize the angle of attack 
flow nonuniformity.. Blade surface fluctuations are seen to be significantly reduced 
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by this streamlining effect. An examination of the overall unsteadiness induced on 
the propfan blades is illustrated by the surface static pressure unsteady envelope 
plot given in Fig. 5.27 corresponding to 90% span. The low level of unsteadiness is 
evident from the thin shaded regions defining the unsteady envelope. The i nfl uence 
of the passage shock clearly dominates this flow and is clearly displayed in the jump 
in static pressure on the pressure surface at roughly 50% chord. An instantaneous 
propfan surface static pressure contour plot is given for the fine grid calculation in 
Fig. 5.28. The asymmetric loading on the spinner and duct leading edges resulting 
from the angle of attack are clearly evident. 

An interesting comparison of unsteady blade loading between an unducted prop- 
fan and a ducted propfan is illustrated in Fig. 5.29. Single blade power coefficient 
rotational history data from the unducted SR7 calculation described in the previous 
section is compared with predicted data from a second calculation for the ducted SR7 
geometry to illustrate the influence of the duct on the time-dependent blade loading. 
In this second ducted propfan calculation, the blade setting angle was increased 3 
degrees to unchoke the flow in the blade passage, and provide a more meaningful 
comparison of results. The overall power coefficient for each case was roughly equiv- 
alent. The solutions indicate that the peak time- dependent loading for the ducted 
propfan is roughly 30% of that predicted for an eqivalently-powered unducted prop- 
fan. This observation suggests that the addition of the propfan duct could provide 
significant reductions in unsteady blade stress levels and propfan generated noise 
levels resulting from angle of attack over comparable unducted propulsion systems. 
Further investigation is required to verify this result over a wider range of geometries 
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and operating conditions. 


5.4 NASA 1.15 Pressure Ratio Fan 

The final calculations to be presented were performed for a 1.15 pressure ratio fan 
stage originally tested by NASA [28]-[31] and utilized extensively under this contract 
for analysis [10]. A description of the geometry and design parameters for the NASA 
1.15 pressure ratio fan is given in Fig. 5.30. This fan is representative of a 25:1 bypass 
ratio turbofan engine fan stage, and therefore closely approximates the ducted propfan 
concept propulsion system. During the present study, only viscous flow calculations 
were performed for this geometry, as a detailed presentation of steady state inviscid 
flow calculations was given in a previous report [10]. For this series of calculations, 
it was not possible to exactly match the experimentally measured mass flow through 
the fan due to the expense of iterating on the three-dimensional solution; therefore, 
predictions are based on the flight Mach number and estimated fan rotational speed 
alone. 

Preliminary steady flow calculations were performed to permit a comparison of 
predicted results with the high speed experimental data published in Ref. [28]. A 
meridional view of the geometry and steady flow viscous flow grid system are given 
in Fig. 5.31. Steady flow calculations were performed for flight Mach numbers of 
0.75 and 0.85. The fan rotational speed was based on advance ratios of 2.86 and 
3.22 for the 0.75 and 0.85 Mach number cases, respectively. A comparison of the 
viscous predicted and experimental cowl leading edge surface static pressure ratios is 
given in Fig. 5.33 for a Mach number of 0.75 and an advance ratio of 2.86. Inviscid 
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Figure 5.25: Full rotor grid system for ducted SR7 propfan geometry 
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Figure 5.26: Instantaneous blade passage static pressure contours for ducted SR7 

propfan geometry (M=0.8). 
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Figure 5.27: Comparison of time-average and unsteady blade section static pressure 

distributions for ducted SR7 propfan (90% span). 
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Figure 5.28: Predicted surface static pressure contours for ducted SR7 propfan ge- 

ometry at angle of attack (M=0.8). 
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Figure 5.29: Predicted single blade power coefficient rotational histories for ducted 

and unducted 8-bladed SR7 propfans at angle of attack (M=0.8). 
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results from a previous task [10] of this contract are also presented on Fig. 5.33 for 
reference. The agreement between experiment and calculation appears to be very 
good for this case. A similar comparison of results is given for a Mach number of 
0.85 and an advance ratio of 3.22 in Fig. 5.34. Again, the viscous results demonstrate 
good agreement with the experimental data in the high gradient leading edge region. 
The deviation observed between the inviscid results and the experimental data near 
the outer downstream edge of the cowl static pressure distribution was previously 
believed to be due either to a shock-induced flow separation or to a wind tunnel 
sidewall interference effect. An examination of the predicted viscous flow velocity 
vectors in this region given in Fig. 5.32 indicates that the flowfield is is not separated 
downstream of the shock, and that the deviation between prediction and experiment 
is likely due to the wind tunnel sidewall aerodynamic interference. 

The final results to be presented are based on a time-dependent viscous flow 
calculation for the NASA 1.15 pressure ratio fan. The flight Mach number was chosen 
to be 0.20, and the advance ratio was 1.12. The angle of attack was 40 degrees. These 
values were selected to permit a comparison with the low-speed experimental data 
for non-zero angles of attack presented in Ref. [29] for the NASA 1.15 pressure ratio 
fan. The actual wind tunnel velocity for this data was reported to be 144 ft/s, which 
corresponds to a Mach number of approximately 0.12. Unmodified compressible flow 
time-marching schemes typically do not perform well for Mach numbers under 0.2, 
and therfore the calculation freestream Mach number was rather arbitrarily increased 
to avoid low Mach number problems in the numerical solution. The results should 
therefore be interpreted with this discrepency in mind. 
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Figure 5.30: NASA 1.15 pressure ratio fan stage geometry (dimensions in cm) 
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Figure 5.31: 


NASA 1.15 pressure ratio fan viscous flow grid system 
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Figure 5.32: Predicted cowl surface velocity vectors in the vicinity of the outer 

surface shock for the NASA 1.15 pressure ratio fan stage geometry 
(M=0.85) 
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The full rotor numerical grid is illustrated in Fig. 5.35 and utilized 488,320 points. 
Relatively large grid spacings were utilized near solid surfaces to reduce the compu- 
tational cost associated with the time-dependent calculation. The experimental case 
corresponding to this calculation was run at 120% design rotor speed, which equates 
to 10,992 rpm. The low flight velocity, high rotational speed, and large angle of attack 
all contribute to the complexity of this test case. 

During the course of this test calculation, some indications of code instability 
due to small radius cells along the centerline upstream of the spinner leading edge 
were observed. These instabilities occasionally led to divergence of the solution, and 
should be considered an area requiring future study. Unfortunately, due to contract 
time and CPU time limitations, it was not possible to run this case to a time-periodic 
solution; however, the results are presented here for completeness. 

An illustration of the instantaneous surface static pressure contours for the time- 
dependent solution is given in Fig. 5.36. The asymmetric loading on the spinner and 
cowl resulting from the highly angled freestream are clearly depicted in the color con- 
tours. A comparison of the time-averaged and unsteady blade surface static pressure 
ratio distribution envelopes are given in Figs. 5.37- 5.39 for 10%, 50% and 90% span, 
on the suction and pressure sides of the blade, respectively. The unsteady loading is 
seen to be largest at the hub and the tip of the blade, presumably due to the flow 
directing influences resulting from the proximity of the hub and cowl surfaces. 

One important aspect of ducted propfan performance is the pressure recovery 
of the inlet duct at angle of attack. This aspect of the calculation is illustrated in 
the comparison of experimental and calculated fan face total pressure ratio contours 
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Figure 5.33: Comparison of viscous and inviscid predicted and experimented cowl 
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Figure 5.35: Full rotor grid system for NASA 1.15 pressure ratio fan 
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Figure 5.36: Predicted viscous instantaneous static pressure contours for NASA 1.15 

pressure ratio fan at angle of attack. (M=0.20) 
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Figure 5.37: Comparison of viscous predicted time-averaged and unsteady envelope 

blade surface static pressure ratio distributions for NASA 1.15 pressure 
ratio fan at angle of attack (10% span, M=0.20) 
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Figure 5.38: Comparison of viscous predicted time-averaged and unsteady envelope 

blade surface static pressure ratio distributions for NASA 1.15 pressure 
ratio fan at angle of attack (50% span, M=0.20) 
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Figure 5.39: Comparison of viscous predicted time-averaged and unsteady envelope 

blade surface static pressure ratio distributions for NASA 1.15 pressure 
ratio fan at angle of attack (90% span, M=0.20) 
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given in Fig. 5.40. The experimental data for this case indicated a significant in- 
crease in inlet total pressure loss for incidence angles greater than 30 degrees. This 
deficiency was attributed to flow separation of the cowl internal flow, as indicated by 
the large decrease in total pressure on the windward side of the fan face total pres- 
sure contours in Fig. 5.40. A comparison of predicted and experimental cowl surface 
static/inlet total pressure distributions along the windward face of the cowl is given 
in Fig. 5.41. The flow separation previously discussed causes the pressure distribution 
to become relatively flat downstream of the highlight. This feature is evident in both 
experimental and predicted pressure distributions. 

The results presented here illustrate the details available from predictions ob- 
tained with the ADPAC analysis code, and the usefulness in assessing off-design 
performance issues associated with ducted propfan propulsion systems operating at 
angle of attack. 
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Figure 5.41: Comparison of experimental and instantaneous viscous predicted cowl 

windward side internal static/inlet total pressure ratio distribution for 
NASA 1.15 pressure ratio fan (M=0.2, angle of attack = 40 degrees) 
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0. CONCLUSIONS 


A time-dependent, three-dimensional Euler/Navier-Stokes aerodynamic analysis 
and grid generation scheme has been developed for the numerical analysis of both 
ducted and unducted propfan flowfields at angle of attack. The underlying multi-block 
discretization scheme utilizes a single H-type grid per blade passage for unducted 
propfans, and a coupled system of five grid blocks utilizing an embedded C-grid 
about the cowl per blade passage for ducted propfans. Aerodynamic predictions 
were verified through comparisons with steady state and time-dependent experimental 
results for an advanced unducted propfan design, and a ducted 1.15 pressure ratio 
fan. Time-dependent calculations for single rotation propfans at angle of attack 
have demonstrated good agreement with experimental data and other predictions. 
The capability of accurately simulating the time-dependent aerodynamics about a 
complete ducted propfan at angle of attack has been demonstrated. 

Several comments are in order concerning the various numerical techniques ap- 
plied in this study. It is apparent that the simple boundary layer dissipation operator 
and algebraic turbulence model are not well suited for the complex vortical flows 
encountered in modern propfan blade designs. The time-accurate implicit residual 
smoothing algorithm can decidedly influence the nature of a solution when a large 
CFL number (and hence, excessive damping) is employed. The accuracy of the anal- 
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ysis can be swayed by additional factors, including the unknown deflected shape of the 
propfan blade, errors introduced through poor grid resolution, turbulence modeling, 
and artificial dissipation. In spite of the known algorithmic deficiencies, the analy- 
sis has successfully predicted the time- dependent flow about ducted and unducted 
propfans at angle of attack, and has demonstrated good agreement with available 
experimental data. 
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